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SUMMARY 


This  is  the  second  year  report  which  is  a  part  of  a  three- 
year  study  on  compliant  robotic  structures.  Such  structures  are 
constructed  of  continuously  flexible  elastomeric  tubes  that 
extend,  bend  or  twist  when  pressurized.  The  motion  of  each  tube 
element  under  pressure  depends  on  its  directional  stiffness, 
achieved  through  the  orientation  of  wall  corrugations  and 
reinforcement.  Tube  elements  placed  in  series  or  parallel  are 
being  designed  as  robotic  fingers  and  arms  that  are  fast-acting 
and  have  potential  payload  to  self  weight  ratios  as  high  as  10/1 
for  laboratory-scale  models  and  up  to  3/1  for  full-scale 
prototype  arms. 

The  first  two  years  of  this  study  have  involved  the 
formulation  and  solution  of  mathematical  models  for  the 
mechanical  behavior  of  single  and  multiple  elastomeric 
structures.  The  formulation  of  these  models  has  been  tempered  by 
the  knowledge  of  the  muscle  morphology  and  function  in  animal 
hydrostats.  In  this  regard.  Chapter  I  describes  fourteen  soft 
animal  parts  (worm  bodies,  cephalopod  appendages,  and  vertebrate 
tongues)  and  the  mechanical  principles  operating  in  these 
structures  as  they  move  during  muscle  contractions.  Chapter  II 
presents  a  general,  nonlinear  mathematical  model  of  orthotrophic 
cylindrical  shells  that  undergo  large  rotations  when  pressurized. 
The  numerical  solutions  may  be  directly  employed  in  the  design  of 
rugged,  lightweight  actuators  to  achieve  wrist  action. 


Chapter  III  deals  with  three  topics :  the  detailed  elastic 
analysis  of  the  load-deformation  behavior  of  bellows;  the  use  of 
this  general  analysis  in  the  design  of  strings  of  bending 
elements  suitable  for  use  as  robotic  arms;  and  a  general  analysis 
of  the  elastica , or  the  finite  deformation  patterns  achieved  in 
cantilevered  element  strings  (or  arms)  under  internal  pressure 
and  end  forces  due  to  its  payload.  This  report  ends  with  a  brief 
description  of  a  novel  compliant,  robotic  arm  based  on  satellite 
bellows  elements,  analogous  to  the  longitudinal  muscles  in  the 
appendages  of  animal  hydrostats. 
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Chapter  I 

MUSCLE  MORPHOLOGY  AND  FUNCTION  IN 
ANIMAL  HYDROSTATS 

Lisa  J.  Croner 

INTRODUCTION 


This  chapter  is  about  muscle,  in  particular,  it  presents  an  overview 
of  a  zoologist's  understanding  of  how  muscles  work  in  a  variety  of 
animals'  soft  parts,  in  the  hope  that  animal  structure  can  serve  as  a 
source  of  ideas  for  the  design  of  flexible  robotic  manipulators.  I  have 
been  asked  by  the  engineers  Involved  In  the  design  of  flexible  manipulators 
to  consider  soft  animal  parts  as  muscular  beams  which  undergo  various 
contortions  when  selected  muscles  contract  To  this  end,  I  devote  the 
first  two  parts  of  the  chapter  to  a  description  of  representative 
cross-sections  of  fourteen  soft,  roughly  cylindrical  animal  parts,  and  to  a 
discussion  of  the  simple  mechanical  principles  operating  in  these 
structures  when  the  muscles  In  them  shorten. 

in  the  third  part  of  the  chapter  l  present  a  brief  discussion  of 
remaining  questions  which  are  likely  to  be  relevant  to  animal-inspired 
flexible  robot  design.  I  present  these  questions  here  to  suggest  subjects 
ripe  for  analysis  in  future  Incarnations  of  this  project. 
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MECHANICS  OF  MUSCULAR  HYDROSTATS 

General  Principles 

Muscle  is  familiar  to  us  as  the  tissue  whose  activity  results  in 
animal  motion.  In  simple  terms,  a  muscle  cell  is  an  elongated  body  which 
contains  a  constant  volume  of  cytoplasm  and  which  has  the  ability  to 
become  shorter  and  fatter  forcefully  when  stimulated  electrically.  This 
shape  change  of  the  individual  muscle  cell  is  converted  into  a  shape  change 
of  a  body  part  when  many  muscle  cells  pull  and  push  on  the  tissues 
surrounding  them.  The  stiffness,  elasticity,  and  other  material  properties 
of  these  tissues  are  thus  extremely  important  in  determining  the  body 
changes  resulting  from  a  muscle  contraction. 

"Muscular  hydrostats"  are  animal  structures  in  which  cell  cytoplasm 
and/or  the  constant  volume  of  soft  tissue  act(s)  as  a  hydrostatic  skeleton 
to  transmit  the  work  of  contracting  muscle  (Kier,  1983).  Worm  bodies, 
cephalopod  appendages,  and  vertebrate  tongues  are  examples  of  muscular 
hydrostats.  The  tissues  comprising  them  are  deformable  and  make  up  a 
constant  volume  of  material.  The  contraction  of  muscles  surrounding  and 
embedded  in  such  tissue  causes  changes  in  the  dimensions  of  the  animal 
structure;  because  the  volume  is  constant,  this  results  in  changes  in  other 
dimensions  of  the  structure.  This  idea  is  now  the  basis  for  the  standard 
analysis  of  movement  in  soft-tissue,  constant  volume  structures 
(Chapman,  1958;  Chapman,  1975;  Clark,  1964;  Kier,  1983;  Kier  and  Smith, 
1985;  Wainwright,1982). 

Two  generalizations  can  be  made  about  the  muscular  hydrostats  I 
will  describe.  These  generalizations  facilitate  the  development  of  a 
model  which  can  be  used  to  explain  the  movements  of  these  structures. 
First,  they  are  all  either  roughly  cylindrical  in  shape  or  consist  of  two 
cylinders  separated  by  a  membrane.  Second,  the  muscles  in  these 
structures  are  found  oriented  in  any  of  three  directions  (Figure  1 ). 

Muscles  may  be  oriented  parallel  to  the  long  axis  of  the  cylindrical 
structure;  these  are  "longitudinal''  muscles.  Muscles  may  be  oriented 
perpendicular  to  the  long  axis  of  the  structure  in  "circular",  "transverse" 
or  "radial"  arrays.  Muscles  may  wrap  around  the  periphery  of  the 
structure  at  an  oblique  angle  to  the  long  axis,  these  are  "oblique"  muscles 
(Kier  and  Smith,  1985). 

Muscular  hydrostats  are  capable  of  making  any  of  four  basic 
movements.  They  may  extend,  shorten,  bend,  and/or  twist.  The 
arrangement  of  muscles  within  the  structures  determines  the  mechanical 
principles  utilized  in  producing  these  movements.  These  principles  have 


been  described  by  various  researchers.  Most  recently,  Kier  and  Smith 
0985)  have  outlined  the  mechanics  of  muscular  hydrostats,  and  the 
following  discussion  reflects  in  large  part  their  views 

Elongation 

Elongation  of  a  soft,  constant  volume  cylinder  occurs  when  its 
cross-sectional  area  decreases.  A  decrease  in  cross-sectional  area  is 
produced  by  the  contraction  of  muscles  perpendicular  to  the  long  axis  of 
the  cylinder.  Hence,  one  would  expect  to  find  circular,  radial,  or 
transverse  muscles  in  soft,  cylindrical  animal  parts  which  are  known  to 
elongate. 

The  contraction  of  oblique  muscles  oriented  at  an  angle  0  greater 
than  54° 44'  to  the  long  axis  of  the  cylinder  will  also  cause  elongation. 
Figure  2  demonstrates  why  this  is  so.  A  constant  volume  cylinder  wrapped 
for  one  turn  by  a  single  oblique  fiber  is  slit  longitudinally  to  display  the 
surface  as  a  flat  sheet,  as  shown  in  Figure  2B.  The  length,  radius,  and 
volume  of  a  right  circular  cylinder  in  relation  to  0  are  respectively 

L  =  D  cose 

r  =  D  sin0/2n 

V  =  nr2L 

where  D  Is  the  length  of  the  oblique  fiber.  Substituting  for  r  and  L,  the 
volume  Is 

V  =  D3sfn20  cos0/4tj. 


The  oblique  fiber  length  is  then 


D  -  (4nV/sin20  cos0)1/3, 


and 


D  /  V '  ^  =  (4n/sin20cos0) '  ^ 


For  a  constant  volume,  0/ V 1  varies  directly  with  D  and  thus  serves  as  a 

"fiber  length  parameter”. 

The  fiber  length  parameter  can  now  be  plotted  against  fiber  angle  0, 
such  a  plot  is  shown  in  Figure  2C.  The  fiber  is  shortest  when  the  fiber 
angle  is  5 4044',  and  increases  in  length  as  the  angle  approaches  0°  or  90° 
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The  shortening  of  a  fiber  oriented  at  less  than  54° 44'  will  increase  the 
fiber  angle,  and  will  thus  result  in  a  decrease  in  the  length  of  the  cylinder. 
The  shortening  of  a  fiber  oriented  at  greater  than  54° 44'  will  decrease  the 
fiber  angle,  thus  resulting  in  an  increase  in  the  length  of  the  cylinder. 
Therefore,  obllgue  musculature  can  contribute  to  elongation  and  shortening 
of  cylindrical  muscular  hydrostats. 

The  amount  and  speed  of  elongation  produced  by  a  given  decrease  in 
diameter  are  affected  by  the  ratio  of  the  length  to  the  width  of  the 
cylinder.  This  is  demonstrated  as  follows.  Consider  a  right  circular 
cylinder  of  constant  volume, 

V  =  TTd2L/4  =  TTd02L0/4 

where  d  and  L  are  the  diameter  and  length  respectively,  and  the  subscript 
zero  denotes  initial  dimensions.  The  diameter-length  relationship  is 
therefore 


<J/d0  =  <L0/L),/2 

This  relationship  is  plotted  in  Figure  3.  Suppose  that  the  initial  state, 
where  d/dQ  =  L/L0  =  I ,  is  designated  as  position  A  on  the  plot,  and  that 

initially  the  diameter  is  one-half  the  length,  or  d0  »  L0/2,  as  shown  by 

cylinder  A.  Then,  if  the  diameter  is  increased  by  a  factor  of  two  (d  -  2d0), 

the  length  must  decrease  by  a  factor  of  four  (L  =  L0/4),  as  shown  by 

cylinder  B.  However,  if  the  diameter  is  decreased  by  a  factor  of  two,  the 
length  increases  fourfold,  as  shown  by  cylinder  C. 

it  is  clear  from  Figure  3  that  a  given  decrease  in  diameter  produced 
by  action  of  radial,  circular,  or  transverse  muscles  causes  a  greater 
change  in  length  In  a  cylinder  with  a  relatively  high  length/diameter  ratio 
than  In  one  with  a  lower  ratio.  The  speed  of  elongation  Is  similarly 
greater  in  cylinders  with  high  length/diameter  ratios.  Therefore,  it  takes 
a  smaller  amount  of  radial,  transverse,  or  circular  muscle  contraction  to 
elongate  a  relatively  long  and  thin  cylindrical  muscular  hydrostat  than  it 
does  to  elongate  a  shorter  and  fatter  one.  One  would  expect,  then,  soft, 
cylindrical  animal  parts  which  elongate  a  great  deal  to  be  long  and  thin. 
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Shortening 

Shortening  is  produced  by  the  contraction  of  longitudinal  muscles  or 
of  oblique  muscles  oriented  at  less  than  54°44'  to  the  long  axis  (Fig.  2). 
The  amount  of  extension  determines  the  amount  of  longitudinal  muscle 
contraction  necessary  to  shorten  a  muscular  hydrostat.  The  longitudinal 
muscles  of  an  extensible  structure  may  have  to  operate  over  a  range  of 
300%  extension  and  contraction.  Obliquely  striated  muscle  cells  are 
thought  to  be  capable  of  operating  over  this  range  (Lanzavecchia  and 
Arcidiacono,  1981). 

The  muscle  cells  of  some  long,  thin,  extensible  structures, 
however,  are  not  obliquely  striated,  and  cannot  function  beyond 
elongations  of  roughly  40%  (Rice,  1973  ).  The  contribution  of  these 
muscles  to  shortening  an  extensible  structure  can  be  increased  if  the 
longitudinal  muscles  are  longer  than  the  structure  itself — in  other  words, 
if  the  longitudinal  muscles  running  through  an  appendage  insert  more 
deeply  in  the  body  of  the  animal  than  does  the  appendage  itself.  For 
example,  a  longitudinal  muscle  1.4  times  the  length  of  the  structure  it 
shortens  will  undergo  a  40%  elongation  during  an  80%  elongation  of  the 
structure. 


Bending 

Bending  of  a  cylindrical  muscular  hydrostat  is  produced  by  the 
simultaneous  contraction  of  the  longitudinal  muscles  on  one  side  of  the 
cylinder  and  the  circular,  radial,  or  transverse  muscles  perpendicular  to 
the  long  axis.  The  material  properties  of  tissues  composing  the  muscular 
hydrostats  1  will  discuss  are  such  that  contraction  of  the  longitudinal 
muscles  alone  would  result  in  a  shortening  and  widening  of  the  cylinder. 
Contraction  of  circular,  radial,  or  transverse  muscles  prevents  the 
diameter  of  the  cylinder  from  increasing,  and  thus  prevents  shortening 
from  occurring.  In  this  case,  contraction  of  longitudinal  muscles  can  only 
shorten  one  side  of  the  cylinder,  and  thus  the  cylinder  bends.  If 
longitudinal  muscles  are  placed  along  the  periphery  of  the  cylinder,  the 
bending  moment  exerted  by  a  given  contraction  of  these  muscles  is  greater 
than  it  would  be  if  they  were  placed  closer  to  the  center  of  the  cylinder. 
Animals  that  have  soft  cylindrical  parts  with  both  longitudinal  muscles 
and  circular,  radial,  or  transverse  muscles  and  with  the  longitudinal 
muscles  located  just  beneath  the  skin  of  the  structure  would  probably, 
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therefore,  include  forceful  bending  as  part  of  the  animals’  repertoire  of 
motions. 

Longitudinal  muscles  located  close  to  the  axis  or  center  of  the 
cylindrical  structure  can  also  contribute  to  bending.  For  a  given  amount  of 
contraction,  they  will  do  so  with  a  smaller  mechanical  advantage  than 
longitudinal  muscles  located  more  peripherally.  They  will  also,  however, 
do  so  more  quickly,  and  thus  with  a  greater  "speed  advantage".  One  would 
expect  that  longitudinal  muscles  in  a  structure  that  bends  very  quickly  and 
without  need  of  great  force  would  be  found  near  the  central  axis  of  the 
structure. 

Twisting 

Twisting  is  produced  by  the  contraction  of  obliquely  oriented 
muscle  fibers.  The  direction  of  twisting  depends  on  the  handedness  of  the 
oblique  muscles.  If  right-  and  left-handed  oblique  muscles  contract 
independently,  both  clockwise  and  counter-clockwise  twisting  is  possible. 
If  they  contract  simultaneously,  the  structure  will  resist  twisting  in  both 
directions. 

As  described  above,  oblique  muscle  contraction  may  also  cause 
elongation  and  shortening  of  a  structure,  depending  on  the  longitudinal 
angle  of  these  muscles.  At  54° 44',  these  muscles  are  at  their  minimum 
length;  contraction  of  muscles  at  this  angle  will  not  cause  shortening  or 
elongation  of  the  cylinder,  but  only  twisting. 

If  oblique  muscles  are  placed  along  the  periphery  of  the  cylinder, 
the  moment  through  which  the  torque  is  applied  for  a  given  contraction  is 
greater  than  if  the  oblique  muscles  are  placed  more  centrally.  Therefore, 
one  would  expect  oblique  muscles  used  for  twisting  to  be  located  close  to 
the  skin  of  the  structure. 


WORM  BODIES,  CEPHALOPOD  APPENDAGES,  AND  VERTEBRATE 
TONGUES 
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The  discussion  so  far  has  focussed  in  a  general  way  on  the 
mechanical  principles  operating  in  cylindrical  muscular  hydrostats.  In 
this  section  I  will  consider  how  these  principles  operate  in  fourteen 
flexible  animal  parts.  These  structures  include  some  of  the  appendages  of 
cephalopods  (the  squid  arm,  the  squid  tentacle,  the  chambered  Nautilus 
cirrus,  and  the  OctopusunN),  some  vertebrate  tongues  (those  of 
pangolins,  lizards,  and  cats),  and  the  bodies  of  some  worms  (horsehair 
worms,  leeches,  the  nematode  Ascaris,  earthworms,  two  genera  of 
polychaetes,  and  the  onychophoran  Peripatudi. 

The  muscle  organization  in  each  of  these  structures  takes  one  of 
three  forms.  The  structures  contain  1  longitudinal  muscles  only,  2) 
longitudinal  muscles  and  transverse,  radial,  or  circular  muscles,  or  3) 
longitudinal  muscles;  transverse,  radial,  or  circular  muscles,  and  oblique 
muscles.  This  represents,  in  a  sense,  a  progression  of  muscle 
organization,  in  which  versatility  of  motion  increases  with  the  addition  of 
each  kind  of  muscle.  This  does  not  represent  an  evolutionary  progression, 
but  simply  a  useful  organizing  principle.  The  structures  will  be  discussed 
with  reference  to  this  progression. 

Structures  with  L  ongitudinai  Muscles  Only 

In  the  generalized  muscular  hydrostat  discussed  above,  the 
contraction  of  longitudinal  muscles  can  cause  shortening  or  can  contribute 
to  bending.  However,  the  absence  of  muscles  whose  contraction  could 
exert  transverse  tension  complicates  matters  in  structures  with  only 
longitudinal  muscles. 

Shortening  requires  subsequent  elongation  to  return  contracted 
longitudinal  muscles  to  their  resting  lengths.  Elongation  results  from  a 
decrease  in  diameter  caused  by  the  contraction  of  radial,  circular,  or 
transverse  muscles.  Elongation  can  also  result  from  elasticity  of  tissue 
deformed  by  shortening.  If  these  muscles  are  absent  and  tissue  elasticity 
is  not  significant,  elongation  cannot  occur.  Shortening  would  therefore  be 
a  devastating  motion  for  such  animals  to  make,  and  one  would  expect  them 
to  devise  mechanisms,  either  neural  or  mechanical,  to  prevent  this. 

Bending  requires  that  some  mechanism  prevent  an  increase  in 
diameter  while  longitudinal  muscles  on  one  side  of  the  hydrostat  contract. 
The  prevention  of  an  increase  in  diameter  need  not  be  an  active  muscular 
process,  but  can  result  from  mechanical  properties  of  tissue  surrounding 


11 


the  hydrostat.  By  preventing  an  increase  in  diameter,  tissue  with 
appropriate  properties  would  also  prevent  shortening. 

Thus  one  would  expect  that  cylindrical  muscular  hydrostats  with 
longitudinal  muscles  only  might  be  incapable  of  shortening,  capable  of 
bending,  and  surrounded  by  material  which  would  prevent  increases  in 
diameter. 

Gordius  (Figure  4)  and  Ascaris  (Figure  5)  represent  two  phyla  of 
worms,  Nematomorpha  and  Nematoda  respectively,  whose  members 
function  with  muscles  oriented  only  longitudinally.  Ascaris  moves  by 
thrashing  its  way  through  viscous  media  (Sherman  and  Sherman,  1976). 
Gordiusswms  or  crawls  by  means  of  whiplike  motions  (Barnes,  1980).  A 
thick  cuticle  surrounds  the  bodies  of  these  worms,  and  it  is  thought  that 
this  cuticle  causes  the  muscles  on  the  dorsal  and  ventral  sides  of  the 
worms  to  act  as  antagonists  (Sherman  and  Sherman,  1976).  In  Ascaris, 
the  contraction  of  longitudinal  muscles  on  one  side  of  the  worm  does  not 
cause  shortening  because  the  thick  cuticle  resists  the  increase  in  width 
associated  with  shortening.  The  cuticle  is  composed  of  several  layers  of 
inextensible  fibers,  presumably  collagen,  which  wrap  in  right-  and 
left-handed  helices  around  the  worm.  These  fibers  are  wrapped  at  an 
average  angle  of  75°30'  to  the  longitudinal  axis  (Harris  and  Crofton,  1957). 
As  described  above  (Figure  2),  hydrostatic  pressure  pushing  against  fibers 
wrapped  around  a  cylinder  at  this  angle  will  produce  an  "elongation  force". 
This  force  is  equivalent  to  that  produced  by  the  contraction  of  oblique 
muscles  wrapped  at  this  angle.  There  will  be  no  increase  in  diameter  of 
the  cylinder,  and  the  worm  will  bend,  stretching  the  longitudinal  muscles 
on  the  opposite  side  to  that  of  the  contracting  longitudinal  muscles. 

Eakin  and  Brandenburger  (1974)  investigated  the  fine  structure  of 
gordian  worms  and  found  that  the  cuticle  surrounding  the  body  consists  of 
as  many  as  36  layers  of  helically  wrapped  fibers.  Although  the  angles  of 
these  fibers  have  not  been  measured,  they  are  probably  greater  than 
54° 44’,  since  this  Is  required  to  resist  an  Increase  In  diameter  when 
hydrostatic  pressure  increases  as  a  result  of  longitudinal  muscle 
contraction.  The  sequence  of  events  producing  the  whlpllke  motion  of 
Gordiuses  likely  the  same  as  those  producing  A  scar  is  thrashing. 

Structures  with  Longitudinal  and  Transverse ,  Radial,  or  Circular 
Muscles 

Structures  with  longitudinal  muscles  and  their  antagonists,  the 
circular,  radial,  or  transverse  muscles,  are  capable  of  elongating. 


Figured.  The  Nematomorph  worm  Gordius,  a)  External  view  of  whole 
animal.  Adapted  from  Barnes,  1980.  b)  Cross-section.  Adapted 
from  Eakin  and  Brandenburger,  1 974. 


Figure  5.  The  Nematode  worm  Ascoris.  a)  External  view  of  whole  animal. 
Adapted  from  Riedl ,  1983.  b)  Cross-section.  Adapted  from  Freeman 
and  Bracegirdle,  1971. 


shortening,  and  bending.  The  structures  I  will  discuss  in  relation  to  this 
musculature  are  the  tongues  of  pangolins  (Man  id),  lizards  (  Varan  us),  and 
domestic  cats  (felid),  and  the  bodies  of  leeches  {Hirudo),  earthworms 
(Lumbricud),  and  one  polychaete  ( Nephtys). 

Hirudo,  Lumbricus,  Nephtys are  all  segmented  worms  of  the 
phylum  Annelida.  Their  bodies  are  divided  into  longitudinal  sequences  of 
segments  separated  from  each  other  by  transverse  septa.  This 
arrangement  allows  each  segment  to  act  as  a  little  hydrostat  independent 
of  other  segments. 

Clark  and  Clark  (1960)  investigated  the  musculature  of  Nephtys,  a 
worm  which  burrows  and  swims  by  means  of  undulatory  waves  of  its  body. 
Figure  6  is  a  drawing  of  two  cross-sections  of  Nephtys.  Figure  6A  shows 
the  musculature  seen  within  a  segment,  and  Figure  6B  shows  the 
musculature  seen  in  a  septum.  The  longitudinal  muscles  are  clearly 
antagonized  by  the  dorso-ventra!  transverse  muscles  within  a  segment; 
these  muscles  prevent  an  increase  in  height  of  a  segment  when 
longitudinal  muscles  contract  to  produce  bending.  It  is  less  clear  what 
muscles  or  structures  prevent  an  increase  in  segmental  width,  although 
Clark  and  Clark  found  no  increase  in  width  during  locomotion.  Nephtys 
does  not  shorten  and  elongate  noticeably  (personal  observation).  Certainly 
there  is  not  an  abundance  of  transverse  muscle  to  produce  elongation. 

///rafo  (Figure  7)  and  Lumbricus  (Figure  8)  are  capable  of  bending, 
shortening,  and  elongating  to  incredible  extents.  They  can  also  move  via 
peristalsis.  The  longitudinal  and  circular  muscles  of  their  body  walls  are 
sufficient  to  produce  these  movements,  as  described  above.  Segmentation 
allows  for  the  independent  action  of  individual  body  segments,  and  this 
allows  for  peristalsis. 

Vertebrate  tongues  are  generally  capable  of  shortening,  elongating, 
and  bending.  They  are  often,  though  not  exclusively,  used  to  gather  food 
and  prey,  and  may  sometimes  be  required  to  maneuver  both  quickly  and 
forcefully.  Pangolins  (Figure  9)  and  lizards  (Figure  10)  are  examples  of 
vertebrates  whose  tongues  have  longitudinal  muscles  located  both 
centrally  and  peripherally.  Centrally  located  longitudinal  muscle  may 
function  to  increase  the  speed  advantage  of  shortening,  and  peripherally 
located  longitudinal  muscle  may  function  to  increase  the  mechanical 
advantage  of  shortening. 

Doran  and  Baggett  (1971)  classified  mammalian  tongues  into  two 
types,  intra-oral  and  extra-oral.  Intra-oral  tongues  are  used  to 
manipulate  food  inside  the  mouth,  and  are  capable  of  less  than  50% 
extension.  Extra-oral  tongues  are  used  in  addition  to  gather  food  outside 
of  the  mouth,  and  are  capable  of  greater  than  100%  extension.  Intra-oral 
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Figure  6.  The  polychaete  worm  Nephtys.  a)  External  view  of  whole 
animal.  Adapted  from  Morris,  Abbot,  and  Haderlie,  1 980. 
b)  Cross-section  showing  musculature  of  septum,  c)  Cross-section 
showing  musculature  of  segment,  b  and  c  adapted  from  Clark  and 
Clark,  1960. 
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Figure  7.  The  leech  Hirudo.  a)  External  view  of  whole  enimal.  Adapted 
from  Barnes,  1980.  b)  Cross- sect  Ion  of  oesophaosal  region, 
c)  Cross-section  of  mid-gut  region,  b  and  c  adapted  from  Freeman 
and  Bracegirdle,  1971. 
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Figure  8.  The  Annelid  earthworm  Lumbricus.  a)  Externa!  view  of  whole 
animal.  Adapted  from  Wainwright  et  al ,  1976.  b)  Cross -sect ion  of 
oesophageal  region,  c)  Cross-section  of  intestinal  region,  d)  Detail 
of  longitudinal  muscle  bundles,  b,  c,  and  d  adapted  from  Freeman 
and  Bracegirdle,  1971. 
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Figure  9.  The  African  pangolin  Mortis,  a)  External  view  of  whole  animal. 
From  Kingdon,  1971.  b)  Cross-section  of  rostrel  part  of  tongue, 
c)  Cross-section  of  middle  part  of  tongue,  d)  Cross-section  of 
caudal  part  of  tongue,  b,  c,  and  d  adapted  from  Doran  and  Allbrook, 


tongues  have  relatively  low  length/diameter  ratios,  while  extra-oral 
tongues  have  relatively  high  length/diameter  ratios. 

An  example  of  an  extra-oral  tongue  is  that  of  the  pangolin  (Figure 
9).  Pangolins  are  African  mammals  that  use  their  tongues  to  dig  for  ants 
and  then  to  manipulate  ants  into  their  mouths.  A  pangolin  tongue  has  a 
high  length/diameter  ratio  and  is  very  extensible.  The  longitudinal 
muscles  of  the  tongue  originate  far  back  in  the  body  of  the  pangolin  and 
stretch  into  the  tip  of  the  tongue.  This  may  be  a  mechanism  for  enabling 
the  longitudinal  muscles  to  retract  the  tongue  after  extension,  as 
discussed  in  the  subsection  "  Shortening “  above. 

Although  Doran  and  Baggett's  classification  was  meant  to  apply  to 
mammalian  tongues,  it  can  be  extended  to  include  many  other  vertebrates' 
tongues.  Lizards  with  tongues  very  similar  to  that  of  the  monitor  lizard 
(  Varanud)  (Figure  10)  have  been  shown  to  have  extensible  tongues  with 
high  length/ diameter  ratios  (Kier  and  Smith,  1985;  Smith,  1984). 

An  example  of  a  mammalian  intra-oral  tongue  is  that  of  the 
domestic  cat  (Fe/id)  (Figure  1 1 ).  Cat  tongues  have  a  low  length/diameter 
ratio  and  do  not  extend  far  out  of  the  mouth. 

Structures  with  Longitudinal  Muscles;  Transverse,  Radial,  or 
Circular  Muscles;  and  Oblique  Muscles 

Structures  with  longitudinal  muscles,  their  antagonists,  and  oblique 
muscles  are  the  most  versatile  of  the  muscular  hydrostats.  They  may  be 
capable  of  shortening,  extending,  bending,  and  twisting.  Such  structures 
include  the  segmented  polychaete  worm  Nereis  (Figure  12),  the 
Onychophoran  worm  Peripatus (Figure  13),  the  squid  Uo/igo)  arms 
(Figure  14)  and  tentacles  (Figure  15),  the  chambered  Nauti/us<z\rn 
(Figure  16),  and  the  arms  of  Octopus  (Figure  17). 

The  arrangement  of  muscles  in  these  structures  directly  reflect  the 
mechanical  principles  operating  in  the  generalized  cylindrical  muscular 
hydrostat.  Longitudinal  and  oblique  muscles  are  located  near  the  periphery 
of  the  structures  to  increase  the  mechanical  advantage  during  bending  and 
twisting.  The  cephaiopod  appendages  also  have  longitudinal  musculature 
located  more  centrally,  implying  that  they  are  capable  of  bending  and 
twisting  with  an  increased  velocity  advantage  also. 


Figure  1 1 .  Cross-section  of  the  tongue  of  the  domestic  cat  Felis. 
Adapted  from  Kier  and  Smith,  1 985. 
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Figure  12.  The  polychaete  worm  Nereis,  e)  External  view  of  whole 
animal.  Adapted  from  Morris,  Abbot,  and  Haderlie,  1980. 
b)  Cross-section.  Adapted  from  Freeman  and  Bracegirdle,  1 97 ! 


v,>. 


»*■'« Vi  ■'>  i 


.  .  .  i 


Figure  13.  Cross-section  of  the  onychophoran  worm  Peripetus  Adapted 
from  Freeman  and  Bracegirdle,  1971. 
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Figure  1 4.  The  squid  Lolipo.  e)  External  view  of  whole  animal.  Adapted 

from  Barnes,  1980.  b)  Cut-away  view  of  arm,  showing  musculature 
From  Kler .  1983. 
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Figure  15.  The  squid  Loligo.  a)  External  view  of  whole  animal  with 

tentacles  extended.  Adapted  from  Barnes,  1980.  D)  Cut-away  view 
of  tentacle,  showing  musculature.  From  Kier ,  1 983. 


Figure  16.  The  chambered  Nautilus,  a)  External  view  of  whole  animal. 
Adapted  from  Barnes,  1 980.  b)  Cut-away  view  of  cirrus,  showing 
musculature.  From  Kier,  1983. 
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Figure  17.  Octopus,  a)  External  view  of  whole  animal,  D)  Cut-away  view 
of  arm,  showing  musculature.  FromKier,  1983. 
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CLOSING  THOUGHTS 

In  the  analysis  presented  above,  muscle  has  been  considered  in  the 
conceptually  simplest  way  possib)e~as  tissue  which  shortens  and  thus 
exerts  tension.  Cross-sections  of  fourteen  different  deformable  animal 
bodies,  appendages,  and  tongues  have  been  discussed  in  terms  of  several 
simple  mechanical  principles  that  must  operate  during  motion  of  these 
structures. 

Several  questions  about  the  mechanics  of  such  structures  remain 
unanswered.  Is  there,  for  example,  a  functional  difference  between  radial, 
circular,  and  transverse  muscle?  In  this  chapter,  these  three  types  of 
muscle  have  been  lumped  into  one  category— musculature  which  is 
perpendicular  to  the  longitudinal  axis  of  a  model  cylinder.  It  is  possible, 
however,  that  circular,  radial,  and  transverse  muscles  differ  in  the  amount 
of  work,  perpendicular  to  the  long  axis,  they  produce  for  a  given 
contraction,  and  that  the  directions  in  which  each  of  these  muscles  pull  is 
functionally  important. 

To  answer  this  question,  detailed  anatomical  data  is  needed,  at 
present,  such  data  is  lacking.  For  example,  do  single  cells  in  circular 
muscle  wrap  all  the  way  around  a  cylindrical  muscular  hydrostat  in  a 
circle,  or  do  cells  overlap  so  that  no  one  cell  goes  all  the  way  around 
although  the  muscle  as  a  whole  does?  If  the  former,  then,  since  the 
circumference  varies  directly  with  the  radius  of  a  circle,  a  40% 
contraction  of  circular  muscle  cells  would  result  in  a  40  %  decrease  in  the 
radius  of  the  circle  enclosed  by  the  muscle;  the  area  of  the  circle  would 
vary  with  the  square  of  the  radius.  If  the  latter,  a  40%  contraction  would 
result  in  some  smaller  decrease  in  the  circumference,  the  radius,  and 
hence  the  cross-sectional  area  of  the  circle  enclosed  by  the  muscles.  The 
same  question  arises  when  analysing  the  radial  and  transverse  muscles, 
and  it  is  evident  that  detailed  anatomical  studies  must  be  performed,  and 
generalizations  and  assumptions  must  be  made. 

To  compare  the  radial  and  circular  muscles  of  a  right  circular 
cylindrical  muscular  hydrostat,  assume  that  circular  muscle  cells  wrap 
completely  around  the  circumference,  radial  muscle  cells  travel  from  the 
center  to  the  edge  of  the  circular  cross-section,  and  the  same  number  of 
cells  of  each  type  must  contract  to  produce  a  given  change  in  radius.  Then, 
since  the  circumference  of  a  circle  varies  directly  with  its  radius,  there 
is  no  difference  between  the  decrease  in  cross-sectional  area  produced  by 
a  given  contraction  of  the  circular  or  radial  muscles. 

In  addition,  the  change  in  area  produced  by  a  given  contraction  of 
radial  or  circular  muscles  will  increase  with  the  size  of  the  circle 
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encompassed  by  the  muscles.  Therefore,  circular  muscles  will  produce  a 
larger  change  in  area  if  they  are  placed  at  the  periphery  of  a  cylindrical 
muscular  hydrostat  than  they  will  if  nearer  the  center,  and  radial  muscles 
will  produce  a  larger  change  in  area  if  they  stretch  from  the  center  to  the 
periphery  of  the  hydrostat  than  they  will  if  they  stretch  only  part  way 
across  the  circular  cross-section.  In  the  animal  structures  presented  in 
this  chapter,  circular  muscles  are  always  found  near  the  periphery  of  the 
structures  (Figures  7,  8,  9,  12,  13,  15,  17),  and  radial  muscles  send 
regular  arrays  of  branches  reaching  as  far  to  the  periphery  as  possible 
(Figure  16). 

Transverse  muscle  is  a  bit  more  difficult  to  analyze.  Transverse 
muscle  cells  that  pass  near  the  center  of  the  circular  cross-section  of  a 
cylindrical  muscular  hydrostat  probably  act  similarly  to  radial  muscles. 
However,  analysis  of  transverse  muscle  cells  placed  as  chords  across  the 
cross-section  nearer  the  periphery  requires  knowledge  of  anatomical 
detail  which  is  not  available.  For  example,  the  vertical  and  horizontal 
muscle  bundles  might  act  in  perpendicular  pairs  to  exert  a  diagonal  force 
on  the  skin  if  both  members  of  a  pair  are  connected  to  the  skin.  It  is 
difficult  to  compare  the  change  in  cross-sectional  area  produced  by 
transverse  muscles  with  that  produced  by  radial  and  circular  muscles 
without  detailed  anatomical  information. 

Another  question  of  interest  is  whether  the  shape  and  arrangement 
of  longitudinal  muscle  bundles  is  functionally  significant.  As  shown  in 
Figure  8D,  the  longitudinal  muscles  of  the  earthworm  are  arranged  as 
many  rows  of  oval  bundles  connected  to  connective  tissue  sheets  which 
run  longitudinally  with  the  muscle.  Many  other  animals  have  longitudinal 
muscle  arranged  seemingly  randomly  in  dense  arrays  of  parallel  roughly 
cylindrical  bundles.  The  significance  of  these  arrangements  has  not  been 
addressed  in  any  literature  of  which  I  am  aware. 

How  can  the  information  presented  in  this  chapter  be  used  by 
designers  of  flexible  robotic  manipulators?  To  date,  engineers  involved  in 
this  project  have  designed  inflatable  plastic  tubes  whose  walls  are 
constructed  so  that  the  tubes  bend,  twist,  elongate,  or  shorten  when  they 
are  inflated.  For  example,  an  inextensible  fiber  oriented  longitudinally  in 
the  wall  will  cause  the  tube  to  bend  when  inflated.  Such  tubes  are  not 
constant  in  volume,  and  they  do  not  contain  shortening  elements;  in  these 
ways  they  are  unlike  muscular  hydrostats.  However,  if  one  were  to  film  a 
partially  transparent  muscular  hydrostat  as  an  opaque  muscle  bundle 
contracted,  using  a  zoom  lens  so  that  the  length  of  the  muscle  remained  a 
constant  in  the  field  of  view,  the  film  would  reveal  something  similar  to 
an  inflating  tube  with  an  inextensible  fiber  element.  The  muscle,  filmed 
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so  that  it  appears  with  a  constant  length,  would  play  the  part  of  the  fiber 
element,  and  the  rest  of  the  structure,  growing  larger  as  the  camera 
zooms  to  maintain  the  length  of  the  muscle,  would  play  the  part  of  the 
inflating  tube.  In  this  way,  inextensible  fibers  in  the  walls  of  inflatable 
plastic  tubes  play  a  similar  role  in  the  tubes  to  that  played  by  muscles  in 
muscular  hydrostats.  The  orientation  and  placement  of  each  determine  the 
motion  of  the  structure  when  the  tube  inflates  or  the  muscle  contracts. 

An  insight  to  be  gained  from  animal  structure  is  that  muscles  (or 
inextensible  fibers)  in  particular  places  and  in  particular  orientations  can 
result  in  particular  motions  when  the  structure  (or  tube)  is  activated.  One 
can  conceive  of  fabulous  flexible  tubes  with  fiber  elements  not  only  in  the 
walls,  but  spanning  the  interior  of  the  tubes,  mimicking  the  muscular 
arrays  and  actions  of  animal  structures. 

Muscle  can  serve  as  an  even  richer  source  of  ideas  if  the  variety  of 
muscle  is  considered.  Although  it  is  tempting  to  ignore  this  variety  and 
conduct  analyses  with  a  generalized  view  of  muscle  contraction,  the  truth 
is  that  muscles  are  not  the  same  from  animal  to  animal  and  even  from 
place  to  place  within  the  same  animal.  Muscle  contracts  because  an 
electrical  impulse,  usually  via  a  nerve,  sets  into  motion  a  sequence  of 
events  which  activates  the  contractile  elements  in  muscle.  There  is  an 
infinite  number  of  variations  possible  in  this  chain  of  events.  Different 
patterns  of  nervous  stimulation  can  give  rise  to  muscle  contractions 
which  vary  in  strength  and  frequency.  The  manner  in  which  an  electrical 
impulse  is  conducted  through  a  muscle  cell  can  determine  the  time-course 
of  a  contraction.  The  arrangement  of  contractile  elements  in  a  muscle  cell 
will  determine  the  force  and  amount  of  contraction.  The  presence  or 
absence  of  mitochondria  in  a  muscle  cell  determines  the  source  of  energy 
for  contraction,  and  thus  affects  characteristics  of  contraction.  An 
engineer  attempting  to  understand  animal  muscle  structure  would  be 
well-advised  to  become  familiar  with  these  variations.  Since  it  would 
require  a  small  book  to  treat  muscle  variation  even  barely  adequately,  I 
refer  the  interested  reader  to  McMahon  ( 1 984),  Wilkie  ( 1 976),  and  Hoyle 
(1983).  These  authors  treat  the  subject  in  detail,  and  their  books  may 
serve  as  sources  of  other  relevant  references.  Consideration  of  the 
differences  between  muscles  will  lead  to  a  more  complete  analysis  of 
animal  structures,  and  may  aid  in  the  choice  of  materials  and  method  of 
control  of  robotic  manipulators. 
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Chapter  II 

FINITE  DEFORMATIONS  OF  NONLINEAR, 

ORTHOTROPIC  CYLINDRICAL  SHELLS 

Gary  Orgill  and  James  F.  Wilson 

INTRODUCTION 

The  behavior  of  cylindrical  shells  under  a  variety  of  loading 
combinations  and  boundary  conditions  has  been  the  subject  of  intense  study. 

The  literature  is  replete  with  classical  and  numerical  solutions  involving 
stress,  deformation  and  buckling  mechanisms.  Most  of  the  work  to  date, 
however,  involves  isotropic,  linear,  elastic  materials  where  deformations  are 
assumed  to  be  small.  Wilson  and  Orgill  (1985)  presented  numerical  results  for 
small  deformations  of  an  orthotropic  thin-wall  right  cylinder  made  of  a 
linear,  elastic  material.  Verma  and  Rana  (1983)  investigated  displacements  of 
a  rotating  concrete  cylinder  with  steel  reinforcing  rods  wound  around  the 
cylinder  as  helices.  Reissner  (1970)  calculated  the  general  expressions  for 
stress  and  strain  in  anisotropic,  thin  walled  cylinders  and  Lekhnitskii  (1963) 
solved  several  problems  involving  linear  elastic  orthotropic  solids.  However, 
the  work  involving  finite  deformations  of  cylindrical  shells,  even  for  linear 
elastic  isotropic  materials,  is  sparse.  The  classical  text  of  Green  and 
Adkins  (1970),  discusses  finite  deformation  analysis  for  shells  from  a  general 
viewpoint.  Leonard  (1967)  solved  the  particular  problem  of  large  deformations 
resulting  from  inflating  a  flexible  isotropic  shell  of  revolution.  Presently 
there  are  several  commercial  computer  codes  employing  numerical  techniques  to 
solve  finite  deformation,  cylindrical  shell  problems. 

As  the  use  of  rubber  and  other  polymeric  materials  become  increasingly 
popular  in  structural  design,  an  analysis  is  needed  that  accounts  for  both 
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geometric  nonlinearities  arising  from  finite  aeformations  and  material 
nonlinearities.  In  applications  using  such  materials,  it  may  De  desirable  to 
incorporate  helical  reinforcement  in  the  cylindrical  shell  or  to  model  helical 
corrugations  as  ortnotropy  of  a  continuum.  A  procedure  that  accounts  for  sucn 

material  and  geometric  nonlinearities,  as  well  as  material  orthotropy  along 
parallel  helices  oriented  at  a  constant  angle  to  the  cylinder's  longitudinal 

axis,  is  developed  herein.  See  Fig.  1.  Parametric  studies  show  how  the 
selection  of  the  orthotropic  properties  will  affect  the  deformation  patterns 
and  how  proper  selection  of  such  properties  can  result  in  an  optimal 
deformation  or  load-carrying  capacity.  First,  the  problem  of  axially 
homogeneous  stress  and  strain  is  considered.  These  results  are  then  extended 
to  include  non-axially  homogeneous  stress  and  strain  that  arise  due  to  end 
constraints.  The  assumption  is  made  throughout  this  study  that  the  stress  and 
strain  in  the  circumferential  direction  are  homogeneous  and  that  the  cylinder 
maintains  a  circular  cross-section  as  it  deforms.  Buckling  loads  and 
post-buckling  behavior  are  not  considered. 

The  loadings  on  the  cylinder  are  internal  pressure,  an  end  load  applied 
along  the  longitudinal  axis  of  the  cylinder,  and  a  pure  torque  about  the 
longitudinal  axis  of  the  cylinder,  applied  separately  or  in  combination.  An 
incremental  analysis  technique  is  used  in  which  small  increments  of  load  are 
applied  to  the  cylinder  at  each  step  and  the  incremental  strain  is  computed  by 
assuming  linearity  over  that  small  range.  The  change  in  section  properties 
such  as  wall  thickness  and  radius  are  computed  from  the  strain  increments  at 
eacn  step,  tnus  accounting  for  the  geometric  nonlinearities.  In  addition  the 
material  properties,  assumed  to  be  known  functions  of  strain,  are  updated  as 
well.  After  each  increment  of  load  is  applied,  the  section  properties  and 
material  properties  are  recomputed.  The  procedure  continues  in  this  fashion 
until  the  desired  load  level  is  attained. 
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PART  I -AXIALLY  HOMOGENEOUS  STRESSES  AND  STRAINS 


Incremental  Stresses 

The  stresses  and  strains  are  assumed  to  be  axially  homogeneous,  which  is 
possible  if  both  ends  of  the  cylinder  are  capped  with  flexible  membranes.  For 
increments  of  internal  pressure  Ap,  end  load  AP,  and  torque  AT,  the  incre¬ 
mental  stresses  at  load  step  n  are  derived  from  elementary  theory.  These 
are: 


Ac 


Ap  R. 


60n  t^ 


AP  Ri  AP 

Aozzn  =  IK-  +  TTTpq 


A°rrn  c  0 


AT  _  AT 


(1) 

(2) 

(3) 

(4) 


where  n  *  i  +1*1,  2,  3,... 

Here  2-j  is  the  polar  section  modulus,  and  R,  and  t^  are  respectively  the  mean 
radius  and  thickness  of  the  cylinder  at  load  step  n.  The  radial  stress  Aorrn 
is  vanishingly  small,  consistent  with  the  thin  wall  assumption. 


Constitutive  Relationships  and  Strain  Transformations 

Consider  the  thin-walled  right  circular  cylinder  of  Fig.  1,  with  the 
orthotropy  defined  by  the  constant  helix  angle  e^.  The  initial  condition  i  = 
0  is  the  no  load  condition.  At  load  step  n  *  i  +  1,  the  incremental  strain 
vector  Aen  is  related  linearly  to  the  incremental  stress  vector  Aon  through 
the  4x4  coefficient  matrix  A,,  assuming  sufficiently  small  load  increments. 
That  is: 

Acn  *  Ai  A on  (i  =  n-1  =  0,1,2,...)  (5) 
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and  T  denotes  transpose. 

The  general  classical  constitutive  equations  for  A-j  as  discussed  by 
Lekhnitskii  (1963)  were  reduced  for  the  present  problem  by  Wilson  and  Orgill 
(1985).  The  components  of  A-j  are  as  follows,  where  the  subscripts  denote, 
respectively,  the  row  position,  the  column  position  and  the  load  step. 
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In  Equations  (8),  E^  and  E^  are  the  tangent  moduli  measured  in  simple  tensile 
tests  and  correspond  to  the  slopes  of  the  stress-strain  curves  for  principal 
strains  and  e2 21  in  the  1  and  2  directions,  respectively.  Similarly, 
is  the  tangent  modulus  for  a  simple  shear  test  at  the  principal  shear  strain 
fl2i*  Poisson's  ratio  v  is  taken  as  invarient  with  respect  to  strain  levels 
clli  and  e22i  as  observed  experimentally  for  rubber  by  Goodyear  (1949).  As 
shown  later,  the  angle  of  orthotropy  8^  changes  with  load  increment. 

The  uniaxial  stress-strain  behavior  with  respect  to  the  principle 
directions  of  orthotropy  is  assumed  in  the  following  form: 


°n  ‘  Eo 

E11  +  B  C11  +  C  E11 

(9) 

°22  -  ECI 

e22  +  B'  C22  +  C'  E22 

(10) 

and  the  behavior  in  simple  shear  with  respect  to  these  principal  directions  is 
of  the  same  form,  or 


a12  =  G0  r12  +  Bg  r12  +  Cg  T12 


(ID 


where  the  coefficients  of  the  strain  are  measured  constants.  Tne  tangent 
moduli  with  respect  to  the  strains  at  level  i  in  this  1-2  coordinate  system, 
found  by  differentiating  Equations  (9)-(ll),  are: 
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These  tangent  moduli  are  used  to  compute  the  elements  of  given  by  Equations 
(8)  for  orientation  6-j . 

The  principal  strains  of  Equations  (12) -( 14)  are  related  to  the  strains 


at  load  level  i  in  the  cylindrical  system  through  the  well-known 
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transformation  equations  (Popov,  1981)  given  by: 

2Elli  ■  Eeef  *  E22i  *  (Ee«1  '  E22t>Eos2ei  *  si"2et 

2e 


22i 


ceei  +  ezzi  +  (ezzi  ’  eeei>cos2ei  ’  Yezi  s1n2ei 


Tl2i  *  (€zzi  +  ceei>sin2ei  +  Ye2i  cos2ei 


(15»a) 

(15b) 

(15c) 


Geometric  Non! ineari ties 

If  the  initial  wall  thickness  is  t0,  then  after  the  first  load  increment 
is  applied  there  will  be  a  corresponding  incremental  radial  strain  Aerri#  The 
updated  wall  thickness  after  the  first  load  increment 

ti  *  t0(l  +  Aerri)  (16) 

Similarly,  after  the  second  load  increment,  the  updated  wall  thickness  is: 

t2  -  ti(l  +  btrrZ)  U7) 

Between  two  successive  load  increments,  it  follows  that  the  updated  wall 
thickness  is: 

*i  a  ti-l(l  +  AErri )  (18) 

In  terms  of  the  initial  wall  thickness  t0.  Equation  (18)  becomes: 

i 

ti  *  to  H  (1  +  Aerrk)  =  to  F^i  (19) 

k*l 

The  factor  represents  the  deviation  of  ti  from  t0  after  load  step  i.  The 
closer  Fti  is  to  unity,  the  less  the  solution  is  affected  by  geometric 
nonlinearity. 

A  similar  argument  is  made  for  changes  in  the  mean  radius.  Equation 
(18),  rewritten  in  terms  of  the  mean  radius  and  the  incremental 
circumferential  strain,  is 

«i  =  Ri _i ( 1  +  Aegoi )  (20) 

Likewise,  Hi  can  be  expressed  in  terms  of  R0,  the  mean  radius  of  the  cylinder 
under  zero  load,  or 


(21) 


Ri  *  R0  n  U  +  Ae68k)=  Ro  Fri 
k*l 

Here,  Fpi  represents  the  deviation  of  R-j  from  R0  after  load  step  i. 

Based  on  the  cross-sectional  area  of  a  thin  ring  given  by  Aj  *  2*  k-j  t-j, 
the  updated  expressions  for  the  area  become: 

Ai  *  A-j_i  (1  +  Aceei )  ( 1  +  Aerr1)  (22) 

i 

Ai  *  2nk0t0  n  (1  +  Aeeek ) ( 1  +  Aerrk)  =  2wR0t0Fai  (23) 

k«l 

Likewise,  the  updated  values  of  the  polar  section  modulus,  initially 
given  by  Z0  *  2nR02t0.  are: 


2i  *  Zi-l(l  +  AeQQi )2( 1  +  Aerri)  (24) 

2i  *  2irR02t0  n  (1  +  Aee6k)2(i  +  Aerrk )  =  2irR02t 0Fp-j  (25) 
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The  corresponding  expressions  for  changes  in  length  of  the  cylinder  are: 

£i  *  *i-l(l  +  Aczz-j )  (26) 

i 

£i  *  lo  ^  (1  +  Ltzzi )  s  *0  Fzi  (27) 

k*l 


Computation  of  Strains 

The  total  strains  are  needed  in  order  to  update  values  of  the  material 
properties  E-j,  Ei',  and  Gj  and  to  compute  cylinder  displacements.  Expressions 
for  these  strains  are  now  deduced.  Consider  the  relationship  between  the  wall 
thicknesses  tj  and  t0  given  in  Equation  (16).  The  total  radial  strain  after 
the  first  load  increment  is: 
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and  after  tne  second  load  increment  is: 
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After  i  increments  of  load  it  follows  that 
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The  total  strains  in  the  circumferential  and  longitudinal  directions,  obtained 
in  a  similar  manner,  are  given  by: 

ceei  *  Fn-  -  1  (31) 


ezzi  *  Fzi  "  1 


(32) 


The  total  shearing  strain  Yezi  1S  simply  the  sum  of  the  incremental  shear 
strains,  as  there  are  no  length  changes  involved.  After  i  increments  of  load, 
this  is 


’ezi  *  k; l ! 


(33) 


The  end  rotation  4>-j  of  the  cylinder  may  be  expressed  in  terms  of  shear 
strain.  For  tne  first  and  second  load  increments,  the  respective  rotations 
are 


♦l  *  '  IT  4Tezl 

0 


(34) 


A*2  m  q  ATez2 


(35) 


The  total  rotation  after  two  increments  of  load  is  the  sum  of  Equations  (34) 
and  (35),  or 


£  £  £  (1  +  4ezzl) 

*2  *  q  AY6zl  +  q  AT6z2  *  q  ^AY6z1  +  (1  +  4cQel)  Ay6z2^ 


(36) 


Extending  Equation  (36)  to  i  increments  of  load  leads  to: 


i-1  1  +  Ac 


♦l  '  £  4Tezk  ra*l  iT^Tii 


o  r(i -1)  k*l 


Deformations  for  Single-Turn  Orthotropy 

Consider  the  segment  of  an  unloaded  cylinder  with  radius  R0  and  length 
i0‘.  The  length  is  selected  such  that  a  line  parallel  with  the  principal 
direction  of  orthotropy  makes  exactly  one  turn  about  the  cylinder  as  shown  in 
Fig.  2.  This  condition  is  expressed  as: 

l'Q  *  2irR0tan  eQ  (38) 

After  the  first  load  increment,  Ro  and  eo  chan9e*  The  resulting 
change  in  geometry  is  given  by  Equations  (20),  (26),  (34)  and  (38),  or 

£' 

■  (2*  -  ♦l)Rl  u"  el  ‘  lo  (1  +  iczzl>  *  <2'  + 

It  follows  that  the  length  for  steps  i-1  and  i  are  given  as: 

N  l  *  +  tan  ei-l  (40) 


i!  *  (2*  +  4>i)Ri  tan  e.  *  (2*  +  AYe?i^Ri  tan  6i  1 

Using  Equations  (40)  and  (41),  a  recursion  relationship  to  compute  is 


deduced  as: 


6'  ■  ^  •»*  Vl^ezi’ 


Increment  Selection  and  Converqence 


The  selection  of  the  appropriate  load  increment  will  depend  on  the 
initial  angle  0O,  the  ratio  of  the  elastic  moduli  E , / E 1  and  E-j/G-j,  the 
magnitude  of  the  load,  and  the  type  of  loading.  As  the  load  increment  is 


increased,  the  computational  effort  is  decreased,  but  the  solution  error  is 
increased.  The  solution  error  after  loading  step  i  is  defined  as 


X  Error 


100( | e.-  c*||2 

[TTYTo 


(43) 


where  e-j  is  the  total  strain  vector  with  components  e@ei»  ezzi »  erri  and  T ez i  • 
The  vector  e*  is  tne  exact  total  strain  vector  for  the  corresponding  loading 
level  containing  the  four  strain  components.  The  percent  error  is  the 
Euclidean  norm  of  the  residual  vector  divided  by  the  Euclidean  norm  of  the 
exact  solution.  The  exact  solution  is  obtained  by  choosing  increasingly 
smaller  loading  increments  until  the  solution  converges. 

Typical  results  of  such  studies  are  shown  in  Table  1  which  gives  values 

of  load  increments  that  can  be  used  over  the  range  of  angle  0  <  60  <  90  deg. 

These  load  increments  were  selected  so  that  the  percent  error  as  defined  never 

exceeds  5%  if  tne  values  of  the  total  nondimensional  load  parameters  do  not 

exceed  their  listed  values.  The  maximum  internal  pressure,  applied  torque  and 

longitudinal  load  are  p,  T  and  P,  respectively.  The  material  is  linear  for 

these  cases,  where  the  only  nonzero  constants  of  Equations  (12)-(14)  are  E0, 

E 1  and  G‘ . 
o  o 


PART  II  NON-AX IALLY  HOMOGENEOUS  STRESSES  AND  STRAINS 

Certain  edge  constraints  give  rise  to  nonuniform  stresses  and  strains  in 
the  loaded  orthotropic  cylinder.  To  account  for  such  nonhomogenity  along  the 
lengtn,  the  cylinder  is  divided  into  J-l  segments  each  of  length  Az0  given  by 

Z 

4Zq  =  (44) 

The  incremental  and  total  stresses  and  strains  are  computed  at  each  point 


Table  1.  Load  Increments  Required  for  a  Solution  Error  of  Less  than  5% 
(Linear  Material ) 


Maximum  value  of 

Ei/Ej 

Increment  in 

Load  Parameter 

Ei/Gi 

Load  Parameter 

P  Ro 

*  0.08 

0  0 

10 

3 

0.001 

H 

O 

• 

o 

50 

3 

0.001 

-J—  *  0.07 

E  tR‘ 

0  0  0 

10 

3 

0.005 

*  0.25 

50 

3 

0.002 

e1«0  * 

0  0  0 

10 

3 

0.001 

0.16 


50 


3 


0.001 


j»l,2,...J  as  shown  in  Fig.  1.  The  boundary  conditions  are  incorporated  when 
the  strains  are  integrated  numerically  for  the  deformations.  The 
computational  procedure  follows  closely  that  of  the  axially  homogeneous  case 
of  Part  I.  However,  quantities  such  as  Ae00.j,  Ftl-,  etc.,  are  now^vector 
quantities  of  dimension  J  with  components  Ae00-jj,  Ft.jj,  etc.  The  added  index 
j  denotes  the  length  from  the  coordinate  origin  as  shown  in  Fig.  1. 

Incremental  Stresses 

The  equation  of  equilibrium  for  an  element  of  a  thin-wall  shell  of 
revolution  at  location  j  and  after  loading  step  i,  deduced  by  Timoshenko  and 
Woinowsky-Kreiger  (1959),  is 

o  o  n . 

eejj  wij  _  _i_ 

Rio  "u  *u  <45> 

Here  is  the  meridional  stress  in  the  shell,  p-jj  is  the  radius  of 

curvature  in  the  meridional  direction,  o00-jj  is  the  circumferential  stress  in 
the  shell,  R-jj  is  the  shell  radius,  p^  is  the  total  internal  pressure,  and  t-jj 
is  the  wall  thickness.  If  the  load  increment  is  sufficiently  small.  Equation 
(45)  can  be  rewritten  in  terms  of  the  section  properties  of  the  previous  step. 


eejj  +  wi  j  „  h 
R(i-l)j  P(i-Dj  t(i-l)j 


The  meridional  stress  (Timoshenko  and  Woinkowski -Krei ger,  1959)  is 


(i-l)j  k(i-l)j  (i "1) j 
By  combining  Equations  (46)  and  (47),  the  circumferential  stress  becomes 


Expressions  for  the  incremental  stresses  for  load  step  i-1  follow  from 
Equations  (47)  and  (48),  or 


A„  . . . 


(49) 


Ac 


eeij 


t(i-l)j  2^(i-l)j 


(50) 


Equations  (49)  and  (50)  are  similar  to  the  internal  pressure  components  of 
stress  given  for  the  axially  homogeneous  case  in  Equations  (1)  and  (2).  Note 
that  for  the  homogeneous  stress  case,  P(i-l)j  *  “  and  o^-jj  *  oQQ ij/2. 

To  compute  strains  and  deformations  in  the  cylindrical  coordinate  system, 
it  is  necessary  to  transform  o^ij  to  that  coordinate  system.  Using  the 
stress-transformation  equations  (Popov,  1981)  yields 

cos  2a,. 


1 


c  . .  -  o  . .  r*7r 


sizm 


■] 


(51) 


zzij  Wij  2  2 

where  tne  angle  that  the  wall  of  the  cylinder  makes  with  the 

longitudinal  axis,  is  defined  as  the  gradient  of  the  radial  deformation,  or 


tan  -  ’"rfcljj 


(52) 


In  employing  the  stress  transformation  of  Equation  (51),  it  is  assumed  that 
alpha  is  sufficiently  small  so  that  the  radial  stresses  in  the  cylindrical 
coordinate  system  are  small,  contributing  a  negligible  amount  to  the  cylinder 
deformation. 

While  the  discussion  in  this  section  has  dealt  only  with  stress  and 


deformation  due  to  internal  pressure,  results  can  be  developed  analogously  for 
end  load  and  pure  torque,  as  well  as  for  various  combinations  of  these  three 
types  of  loads. 


Constitutive  Relationships,  Material  and  Geometric  Nonlinearity 

The  incremental  constitutive  law  given  by  Equations  (5 )-(b),  along  with 
the  material  nonlinearity  of  Equations  (12) -( 14)  and  the  transformations  of 
Equations  (15),  are  assumed  to  be  valid  not  only  at  each  load  step  i  but  also 
at  each  length  location  j.  With  this  notation.  Equations  ( 5 ) - ( B )  become 


Ae86nj 

a  1 1  i  j 

a12i  j 

a13i  j 

“ 

a  1 4  i  j 

Ao88nj 

AEzznj 

= 

a21i  j 

a22i  j 

a23i  j 

a24i  j 

Aozznj 

AErrnj 

a  3 1  i  j 

a32i  j 

a33i  j 

a34i  j 

Aorrnj 

AE8znj 

a41i  j 

a42i  j 

a43i  j 

a44i  j 

Aoeznj 

n  =  i  +  1  =  1,2,3,... 
j  *  1,2,3, ...J 


Note  the  components  of  the  coefficient  matrix  of  Equation  (53)  are  still  given 
by  Equations  (8)  where  each  parameter  with  the  subscript  i  now  has  the  added 
subscript  j.  In  Equations  (12)-(15),  the  single  subscript  i  is  also  replaced 
by  ij. 

The  geometric  properties,  except  for  length  change,  are  deauced  from 
Equations  (19),  (21),  (23)  and  (25),  i  being  replaced  by  ij.  The  change  in 
length  along  the  longitudinal  axis  for  a  segment  bounded  by  points  j  and  j+1 
is  given  by 

Azij  =  Azo  +  uzi(j+l)  "  uzij  5  Azo  ^zij 
where  4z0  is  given  by  Equation  (44)  and  where  u2i(j+i)  amd  uzij  are  the  total 
displacements  measured  relative  to  the  initial  points  j  and  j+1  when  the 
cylinder  is  without  load  (i=0). 


Computation  of  Strain,  Rotation  and  Orthotro 


The  total  strains  can  be  computed  either  from  the  incremental  strains  or 
from  the  displacements.  The  radial,  circumferential  and  shear  strain  at  load 
level  i  and  position  j,  expressed  as  Cprij*  eeeij»  and  ^ezij  respectively,  are 
given  explicitly  by  Equations  (30),  (31)  and  (33),  where  ij  replaces  i. 

An  efficient  method  for  calculating  eZZij,  the  longitudinal  strain,  is  by 
the  finite  difference  method  used  by  Utku  (1981).  If  uZij  denotes  a  segment 
displacement,  then 

c  ..  .  "*»♦!)  '  yu-l)  (55) 

The  average  longitudinal  strain  for  the  whole  cylinder  of  initial  length  £0  is 

'czzij  "I"  (uzU  ‘  “zil1  (66) 

where  J-l  is  the  number  of  length  segments. 

The  end  rotation  <p-jj  of  the  cylinder  segment  between  locations  j  and  j+1 
can  be  computed  from  Equation  (37),  after  the  subscript  j  is  added  where 
appropriate;  or  from  tne  expression  involving  the  incremental  ci rcumferential 
displacement  Aue-jj  given  by 

Au.  .  . 

(57> 

The  angle  of  orthotropy  0-j  j  in  each  of  the  J-l  longitudinal  segments  is 
given  by  Equation  (42)  where  each  subscript  i  or  (i-1)  is  replaced  by  ij  or 
(i-l)j,  respectively. 


Computation  of  Displacements  and  Curvature 

The  computation  of  displacements  ur,  uz  and  u@  in  the  case  of  axially 
homogeneous  strain  is  a  straightforward  application  of  the  strain-displacement 
equations  (Wilson  and  Orgill,  1985).  However,  when  the  strains  are  axially 


nonhomoyeneous,  displacements  are  more  difficult  to  compute  and  in  general 
must  be  obtained  numerically.  The  results  that  follow  will  be  based  on  the 
finite  difference  method  discussed  by  Utku  (1981)  and  Carnahan,  et  al  (1964). 

From  the  strain-displacement  equations  given  by  Sokolnikoff  (1966),  the 
incremental  radial  displacement  becomes 

^rij  *  R(i-l)j  Ae66ij  (58) 

The  incremental  displacement  is  then  added  to  the  total  displacement  from  the 
previous  load  step  to  obtain  the  total  radial  displacement  at  step  i  for 
segment  j,  or 

uri j  *  ur(i_i)j  +  Aurij  (59) 

Computation  of  the  incremental  longitudinal  displacement  is  obtained  from 
the  strain-displacement  equations  by  integrating  the  incremental  longitudinal 
strain  given  by 

Auzi(j+1)  -  Auzi ( j-1 )  *  [^(i-ljj  +  Az(i-l)(j-l)3  Aezzij  (60) 

As  j  varies  from  1  to  J  in  Equation  (60),  a  tridiagonal  system  of  linear 
equations  is  found.  These  equations  are  solved  simultaneously  for  the 
incremental  displacements  Auz-jj.  Tne  total  longitudinal  displacement  at  step 
i  is  once  again  computed  by  adding  the  incremental  displacement  to  the  total 
longitudinal  displacement  at  the  previous  step. 

uzij  *  u2(-j_i)j  +  Auz-jj  (61) 

The  computation  of  the  circumferential  displacement  is  obtained  from  the 
strain-displacement  equations  by  integrating  the  incremental  shear  strain. 

The  integration  results  in  expressions  similar  to  Equations  (60)  and  (61)  for 
the  longitudinal  displacement.  The  incremental  circumferential  displacement 
and  total  circumferential  displacement  are,  respectively: 


Auei(j+i)  -  Auei(j-i)  =  Uz(i-i)j  +  Az(i-i)(j-i)]  Ayezij 
u ©i j  *  ue(i-i) j  +  4ueij 


(62) 
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With  these  displacements  it  is  possible  to  compute  the  angle  between 
the  meridional  direction  and  the  longitudinal  axis  of  the  cylinder;  and  p-jj, 
the  radius  of  curvature  in  the  meridional  direction.  The  angle  a^j  is  used  in 
transforming  stress  from  the  meridional  direction  to  the  longitudinal 
direction  as  given  in  Equation  (51).  Equation  (52)  expresses  o.jj  in  terms  of 
the  derivative  of  ur-j j  in  the  z-direction.  In  finite  difference  form,  the 
first  and  second  derivatives  of  this  displacement  are  as  follows. 


~Url (j-1)  *  uri (j+1) 

^(i-DJ  +  ^(i-DU-D 


tan  a. 


(64) 


3u. 


i'Vii  ,  2[uri(j-l)  *  "rt(J+l)  '  2unj  ^'“(i-DIJ-D  ~ 

-2  a «* 2  a -.2  (65) 
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Sz‘  “(t-DU-l)  T 

The  radius  of  curvature  is  computed  from  the  following  familiar  form. 


X 

AZ  #  _•  i  \  /  _•  i\  +  Az  /  ^ 


together  with  Equations  (64)  and  (65). 
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A  comprehensive  treatment  of  the  finite  difference  representation  of 
displacement  boundary  conditions  that  satisfy  various  end  restraints  at  j=l 
and  J  is  given  by  Carnahan,  et  al  (1964).  A  particular  example  is  presented 
below. 

For  suitable  convergence,  it  is  important  to  select  a  sufficient  number 
of  cylindrical  segments  J-1  to  accurately  represent  the  cylinder's  overall 
shape.  As  Equations  (60)  and  (61)  are  integrated  to  yield  the  displacement, 
the  error  in  the  result  will  decrease  as  J  is  increased.  However,  if  J  is  too 
large,  the  accuracy  of  the  derivatives  computed  by  Equations  (64)  and  (65)  may 
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actually  decrease.  It  was  found  that  the  selection  of  the  number  J  is 
tempered  by  the  initial  cylinder  geometry,  the  orthotropy  angle  e0,  the 
material  constants,  the  type  and  magnitude  of  loading,  and  the  computational 
effort. 


PARAMETRIC  STUDIES 

A  Fortran  IV  computer  code  was  written  to  carry  out  the  calculations  for 
the  finite  deformation  behavior  of  both  the  homogeneous  and  nonhomogeneous 
types  of  cylinders.  Selected  parametric  studies  for  both  types  are  now  pre¬ 
sented.  For  cases  where  the  material  is  linear  and  only  one  of  the  three 
loads  is  present  (longitudinal  load  P,  internal  pressure  p,  or  pure  torque  T), 
the  deformations  depend  on  the  following  appropriate  independent  nondimen- 
sional  load  parameters  and  three  geometric  parameters. 
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(67) 


In  the  selected  examples,  the  deformations  were  found  to  be  very  sensitive  to 
tne  independent  parameter  e0,  tne  initial  angle  of  orthotropy.  Thus  60  was 
taken  as  the  abscissa  in  the  presentation  of  the  graphical  results.  The  last 
two  geometric  parameters  of  (67)  affect  deformations  only  in  axial 
nonhomogeneous  problems. 

The  first  type  of  cylinder  has  axially  homogeneous  stresses  and  strains 
and  is  made  of  an  orthotropic,  linear  elastic  material  with  the  following 
properties: 


E-j/E-l '  -  E0/E0'  «  10  ;  Ei/Gi  =  E0/G0  =  3 
B  *  B*  *  Bg  *  C  *  C  *  Cg  8  0  ;  v  =  0.5 
The  results  for  this  cylinder  are  presented  in  Figs.  3-8. 


(68) 


Figures  3,  4  and  5  show  the  behavior  of  the  finite  longitudinal  strain 
ezz  with  the  loading  P,  p,  and  T  respectively.  For  longitudinal  loading,  it 
is  deduced  from  Fig.  3  that  for  0O  *  27  deg  and  60  >  60  deg,  the  curves 
coalesce  and  ezz  is  linear  with  P.  For  0  <  60  <  27  deg,  however,  the  cylinder 
acts  as  a  "soft’1  spring  in  the  axial  direction,  since  as  P  increases,  ezz 
becomes  proportionally  larger.  For  27  <  60  <  60  deg,  however,  the  reverse  is 
true  but  is  less  pronounced,  where  the  cylinder  behaves  as  a  "hard"  spring. 
However,  for  internal  pressure  loading,  Fig.  4,  linear  behavior  is  found  for 
e0  ■  52  deg,  with  soft  spring  behavior  for  smaller  values  and  hard  spring 
behavior  (but  less  pronounced)  for  larger  values  of  60.  Except  near  the 
extreme  values  of  e0,  the  behavior  of  axial  strain  with  pure  torque  is  more 
complex,  as  shown  in  Fig.  5.  For  0O  up  to  3  deg,  the  behavior  is  linear;  and 
from  that  point  to  30  deg  the  cylinder  “winds  together"  where,  as  T  increases, 
changes  in  ezz  are  proportionally  smaller.  However,  for  e0  near  80  deg,  as  T 
increases,  ezz  becomes  proportionally  larger,  and  the  cylinder  winds  together 
more  and  more  easily,  witnin  small  cnanges  for  0-j.  It  is  noted  that  0-j 
decreases  as  the  cylinder  winds  together. 

Figures  6,  7  and  8  show  tne  behavior  of  $,  the  finite  angle  of  rotation 
about  the  longitudinal  axis  of  the  cylinder,  with  loading  P,  p,  and  T 
respecti vely.  The  negative  ordinate  in  Figs.  6  and  7  indicates  that  the 
cylinder  unwinds  with  the  application  of  load.  The  change  from  a  soft  or 
flexible  configuration  to  a  hard  one  at  0O  *  30  deg  is  apparent  in  Fig.  6;  but 
for  pressure  loading  this  change  is  less  distinct  and  occurs  when  e0  is 
between  50  and  60  deg,  as  shown  in  Fig.  7.  The  effect  of  pure  torque  on  $  for 
0O  up  to  about  30  deg  is  clear  from  Fig.  8:  as  T  increases,  $  becomes 
proportionally  larger,  and  the  cylinder  winds  together  and  becomes  shorter. 
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Effect  of  axial  load  and  ortnotropy  on  axial  strain  (axially 
homogeneous  case) 
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Effect  of  axial  load  and  orthotropy  on  end  rotation  (axially 
homooeneous  case) 
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The  second  type  of  cylinder  selected  for  study  has  axially  nonhomogeneous 
stresses  and  strains,  but  has  the  same  material  constants  as  for  the  study 
just  completed.  That  is,  the  material  is  orthotropic,  linear  and  elastic  with 
the  properties  defined  by  Equations  (68).  For  axial  nonhomogeneity,  the 
initial  geometry  (£0,  R0,  t0)  and  number  of  points  J  along  the  length  affect 
the  deformation  patterns.  These  particular  quantities  were  chosen  as 
follows: 

VR0  *  10  ;  V^o s  10  ;  J  •  n  (6s> 

The  boundary  conditions  for  the  second  type  of  cylinder  were  chosen  to 
simulate  a  fully  clamped  condition  at  one  end,  j*l,  and  zero  radial 
displacement  at  the  other  end,  J  *  11.  However,  the  longitudinal  displacement 
and  twisting  angle  were  chosen  to  be  unrestrained  at  J  *  11.  These  conditions 
are: 

u2i  *  uel  *  url  *  urj  *  0  (70) 

The  results  for  this  second  type  of  cylinder  are  presented  in  Figs.  9,  10 
and  11  in  which  the  only  loading  is  internal  pressure.  For  the  range  of  p 
that  overlaps  (0.001  <  p  <■  0.02),  it  is  observed  that  the  curves  of  Fig.  9  for 
the  average  longitudinal  strain  e22  are  identical  to  those  of  Fig.  4  for  tne 
homogeneous  strain  e22.  For  this  range,  then,  the  end  boundary  constraints 
have  a  negligible  effect  on  the  total  longitudinal  deformation.  However,  for 
the  range  of  p  that  overlaps  (0.001  <  p  <  0.02),  the  total  end  rotations  4>  of 
Fig.  10  differ  from  the  counterpart  results  of  Fig.  7  (boundaries  uncon¬ 
strained),  where  these  rotations  become  more  dissimilar  as  ©0  is  increased 
beyond  15  deg.  Compared  to  the  constrained  case,  the  geometric  changes  along 
the  length  of  an  unconstrained  cylinder  are  always  more  pronounced.  The  last 
study.  Fig.  11,  shows  the  nonlinear  variation  in  midlength  radius  R^,  with 
internal  pressure.  Such  calculations  are  of  particular  importance  for  axially 
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nonhomogeneous  deformations  because  they  serve  as  a  way  to  check  the  limits  of 
validity  of  the  present  analysis.  That  is,  since  out-of-plane  stresses, 
strains,  and  deformations  were  ignored,  then  the  angle  o^j  needs  to  be 
sufficiently  small  for  the  results  of  Figs.  9,  10  and  11  to  be  valid,  or  cos 
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CHAPTER  III 


i 

LARGE  DEFLECTIONS  OF  CONTINUOUS  ELASTIC  STRUCTURES 
Murugappan  Palaniappan  and  James  F.  Wilson 

1.  INTRODUCTION 

A  robot  is  a  reprogrammable,  multifunctional  manipulator  designed  to 
perform  a  variety  of  tasks.  Examples  of  robot  functions  are  material 
handling,  spray  painting  and  aligning  of  screws  in  assembly  line  operations 
(Critchlow,  1985).  The  present  state  of  robot  technology  is  that  the  arm 
moves  around  accurately  but  not  efficiently.  Most  of  the  robotic  arms  in  use 
today  consist  of  rigid  elements  connected  by  hinged  or  pivoted  joints.  The 
movement  of  these  arms  is  slow  because  of  their  high  weight  and  inertia,  and 
the  time  needed  to  compute  and  control  the  coordinates  of  their  elements  or 
arm  segments. 

A  different  approach  to  arm  design  is  based  on  the  motion  of  a 
flexible  elephant  trunk  (Mahajan,  1985).  The  flexible  arm  is  made  of 
lightweight,  polymeric  material  and  may  handle  comparatively  heavy  payloads  at 
higher  speeds  than  most  rigid  arm  designs.  These  flexible  arms  are  based  on 
the  action  of  corrugated  tube  elements  (Wilson,  1984-a)  that  bend  when 
pressuri zed. 

Three  distinct  aspects  of  flexible  arm  design  are  analyzed  in  this 
thesis:  the  equivalent  or  reduced  modulus  E'  for  the  corrugated  portion  of  a 
bending  element;  the  single,  elastic  bending  element  treated  as  a  composite, 
reinforced,  elastic  cantilevered  beam;  and  the  whole  arm  deflection  analysis, 
treated  as  bending  elements  in  series.  In  this  analysis,  self  weight  is 


neglected  since  the  payloads  are  expected  to  be  at  least  an  order  of  magnitude 
higher  than  the  arm  weight. 

Section  2  deals  with  the  computation  of  reduced  modulus  E‘  of  a 
rectangular  corrugation  in  which  the  static  load  responses  of  bellows  type 
arms  with  Young's  modulus  E»  using  classical  theory  of  plates  and  shells  are 
derived.  In  his  classic  paper,  Donnell  (1932)  calculated  the  longitudinal 
flexibility  for  several  types  of  corrugations  and  derived  a  reduced  modulus  of 
elasticity  for  these  corrugated  bellows  as  if  they  were  smooth  tubes. 

Calladine  (1974)  used  energy  methods  to  compute  the  flexibility  of  an  axially 
symmetric  elastic  bellows  subjected  to  axial  loading.  Haringx  (1952)  studied 
the  instability  of  rectangular  corrugation  bellows  when  loaded  by  internal 
pressure.  Wilson  (1984-b)  reviewed  several  mathematical  models  and  the 
assumptions  used  to  compute  E'  and  to  compare  calculated  and  experimental 
values  of  E'.  He  also  shows  that  the  mean  bellows  radius  is  the  appropriate 
dimension  needed  to  convert  internal  pressure  to  an  external  load. 

Section  3  deals  with  the  plane,  finite  deflection  analysis  of  a 
cantilever  beam  subjected  to  end  loading.  The  exact  differential  equation  is 
solved  in  terms  of  the  slope  of  the  elastic  curve.  Related  work  is  summarized 
as  follows:  Scott,  Carver  (1953)  presented  an  integral  power  series  solution 
of  the  nonlinear  beam  equation  in  which  the  moment  is  expressed  as  a  function 
of  the  distance  x  from  the  origin.  Lewis,  Monasa  (1982)  presented  a  large 
deflection  analysis  of  thin  cantilever  beams  of  nonlinear  materials  subjected 
to  a  constant  end  moment.  Theocaris,  Panayotounakos  (1982)  solved  the 
nonlinear  differential  equation  of  an  elastic  cantilever  subjected  to  coplanar 
terminal  loading,  taking  into  account  the  influence  of  transverse  shear 
deformations.  Bisshop,  Drucker  (1945)  presented  a  solution  for  large 
deflections  of  a  cantilever  beam  subjected  to  a  tip  transverse  load  in  terms 


of  an  elliptic  integral. 

The  solutions  given  by  ([Scott,  Carver],  [Lewis,  Monasa],  [Theocaris, 
Panayotounakos])  are  in  terms  of  Cartesian  coordinates.  The  solution 
presented  herein  involves  solving  iteratively  a  complete  elliptic  integral  of 
the  first  kind  for  the  maximum  slope  at  the  tip  of  an  elastic  cantilever 
subjected  to  general  plane  loading  at  the  tip.  This  solution  is  a 
generalization  of  the  elastica  in  which  only  longitudinal  loads  were 
considered  (Timoshenko  and  Gere,  1961). 

Section  4  deals  with  the  deflection  of  the  flexible  arm  in  two 
dimensional  space,  where  the  arm  is  composed  of  end-to-end  cantilevers.  Each 
element  of  the  arm  is  analyzed  as  an  elastic  cantilever,  starting  from  the 
tip,  with  the  fixed  end  at  the  tip  of  the  previous  element.  Transformation 
equations  are  used  to  relate  the  tip  of  each  element  to  the  global  coordinate 
system  which  is  fixed  at  the  root  of  the  arm.  Knowing  the  final  destination 
of  the  payload  at  the  tip  of  the  arm,  an  iterative  procedure  is  used  to  vary 
the  pressure  history  in  each  corrugated  tube  element  of  the  arm  so  that  the 
payload  may  reach  its  destination. 

In  Section  5,  conclusions  of  this  study  and  recommendations  for 
further  research  are  presented.  Different  geometrical  arrangement  of  the 
bellows  is  discussed  to  achieve  out-of-plane  motion.  A  systematic  procedure 
is  suggested  to  vary  the  pressure  to  achieve  rapid  convergence  to  the 
destination  coordinate.  A  data  base  approach  is  suggested  to  overcome  the 
numerous  computation  involved  in  the  computation  of  the  tip  angle  a. 

All  computer  codes  listed  in  the  appendices  were  written  by  the 
author.  They  are  used  to  compute  the  reduced  modulus  E',  the  shape  of  the 
elastic  curve  of  the  deformation  behavior  of  a  single  cantilever  and  the 
equilibrium  states  of  a  manipulator  arm. 


2.  ANALYSIS  OF  BELLOWS  ELEMENT 


A.  Mathematical  Model 

Bellows  are  thin  walled  corrugated  tubes  designed  for  high  flexibility 
when  subjected  to  axial  loads,  internal  pressure  or  bending  moments.  The 
present  study  is  of  a  bellows  corrugation  model  consisting  of  cylindrical 
shell  sections  connected  by  annular  plates  and  subjected  to  axial  load  and 
internal  pressure.  These  bellows  are  to  be  used  as  the  structural  section  of 
the  robotic  manipulator.  The  desired  motion  of  the  bellows  will  be  controlled 
by  microcomputer  and  monitored  through  pressure  gages.  Figure  (2.1)  shows  a 
longitudinal  section  through  the  bellows,  which  is  subjected  to  internal 
pressure  p  and  axial  force  F0.  The  shaded  section  is  a  typical  half 
corrugation  used  in  the  analysis. 

Calculations  for  bellows  extension,  compression  and  transverse  bending 
are  made  as  if  the  bellows  were  an  equivalent  cylinder  of  the  same  mean 
radius,  length  and  wall  thickness.  The  equivalent  cylinder  has  a  reduced 
Young's  modulus  E'  which  is  used  to  compute  static  load  responses  of  the 
bellows  with  modulus  E.  Figure  (2.2)  shows  the  geometry  of  a  half  corrugation 
of  the  bellows  together  with  the  coordinate  system  for  the  inner  and  outer 
shells  and  plate.  For  a  bellows  of  n  corrugations  of  mean  cross  sectional 
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Figure  2.2  Geometry  of  half  corrugation  with  local  coordinate  system. 


area  2irRt,  the  axial  extension  from  strength  of  materials  is 

P  4nb  _  P  2nb 
*b  =  E1  2  irRt  "  V  irRt 


(2.1) 


where 

P  =  F0  +  irR^p 

E'  =  modulus  of  equivalent  cylinder 

As  the  corrugation  is  symmetric  transversely,  we  need  to  consider  only 
one-half  of  the  corrugation,  of  length  2b,  to  transform  it  into  an  equivalent 
cylinder,  as  shown  in  Figure  (2.3).  The  axial  extension  of  the  equivalent 
cylinder  is 


b 

irRt 


(2.2) 


or 


First,  the  axial  extension  A  of  the  corrugation  is  predicted  from  the 
theory  of  plates  and  shells.  Equation  (2.2a)  is  then  used  to  compute  the 
reduced  modulus  E‘  of  the  equivalent  cylinder.  Figure  (2.4)  shows  a  free  body 
sketch  of  the  corrugation  together  with  its  slope  and  moment  compatibility  at 
the  junctions.  Figure  (2.5)  shows  the  corrugation  subjected  to  internal 
pressure  p  and  axial  force  F0  and  its  resulting  deflected  shape.  Radial 
displacement  of  the  plate  is  neglected. 

The  assumptions  of  the  mathematical  model  are  summarized  below: 

1)  The  inner  and  outer  cylindrical  shells  restrain  the  radial  displacement 
of  the  annular  plate  at  junctions  1  and  2. 

2)  Right  angles  are  maintained  at  junctions  1  and  2. 

3)  Axial  forces  in  the  cylindrical  shells  have  negligible  effect  on  the 


shell  stiffnesses. 


MODULUS  E 


MODULUS  E 


Figure  2.3  Half  corrugation  of  modulus  E  with  axial  elongation  A  and 
its  equivalent  cylinder  of  modulus  E’  with  same  axial  elon¬ 
gation  A. 


4)  The  cylindrical  shells  contribute  a  negligible  amount  to  the  axial 


deflection  under  internal  pressure  p  and  axi symmetrical  load  F0. 

5)  Plate  bending  accounts  for  all  of  the  axial  deflection  under  loads  p  and 

F0. 

6)  In-plane  (radial)  loads  on  the  plate  have  negligible  effect  on  its 
stiffness  and  axial  deflection. 

7)  Effects  of  rigid  end  plates  on  symmetrical  deformation  of  a  typical 
element  are  negligible. 

8)  Classical  plate  and  shell  theory  govern  the  system  behavior  under  loads  p 
and  F0. 


B.  Plate  Analysis 


The  pressure  loading  p  is  a  non-negative  quantity  and  the  sign  associated  with 
p  is  negative  as  the  direction  of  p  is  opposite  to  the  direction  of  lateral 


deflection  w.  Here,  r  is  the  radial  coordinate  and  Dp  is  the  rigidity  of  the 
plate. 

A  more  convenient  form  of  Equation  (2.3)  is 


(2.4) 


7  77^777^>  =  -f 

P 

The  moment  and  shear  equations  of  a  plate  are 


u  _  n  /d  w  v  dwv 

Mrp  '  *  °p  ^  +  7  dF> 

0  =  -  D  4—  [—  ^—  ( r  $!■)] 

rp  p  dr  r  dr  dr 
Integrating  Equation  (2.4)  directly,  we  obtain 


(2.5) 


(2.6) 


*  ■  -  ♦  C3  (%)2  "  <%>  *  C2  %>'  +  C1W  (%>  *  Co 

where  CQ,  Cj,  C2  and  C3  are  arbitrary  constants  of  integration. 

These  four  constants  are  written  in  terms  of  the  following  edge 


~)  +  C 


(2.7) 


conditions. 


w  =  0  at  r  =  r. 


(2.8) 


u  _  n  /d  W  ,  V  dWv 

7  ^ 


(2.9) 


m  _  pi  /d  w  ,  v  dw \ 

%p  °o(^ 

%  ■  -  °p  Sr  C7  Sr  <r37'] 

2p 


Apply  Equation  (2.8)  to  (2.4). 


(2.10) 


(2.11) 


rl  p 

Co  =  64Dp  "  C2 

Substitute  C0  in  Equation  (2.4). 


'  ■  -  <r4  -  rl4>  A-  +  <?7>  *  C2  ^  -  »  +  Clw  (77’  <2' 


r  x2„_  , r 


(2.12) 


T-)  (2-13) 


Use  the  derivatives  of  w  in  Equation  (2.9)  and  evaluate  the  moment  at 


73 


p  ^*p  c  c  c 

M_  =  0  {-  ■  (3  +  v)  +  (3  +  v)  +  -^2  2(1  +  v)  -  -^2  (1  ■  v) }  (2.14) 


Plp 


1 


1 
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Use  the  derivatives  of  w  in  Equation  (2.10)  and  evaluate  the  moment  at 


r  *  r«. 


or 


M 


r2p 


°P  {- 


r2  P 

TSd7 


(3  +  v)  +  — *-  [3  t  v  +  2  in  (-^-)(1  +  v)] 

r  c  ri 

P1  1 


c  c 

+ -2^2(1  +  v)  --±2  (1  -  v)} 


'1  '2 

Use  Equation  (2.11)  and  evaluate  the  shear  at  r  *  r^, 

,r2P 


'  r 

r2p 


(20“  -  C3  ~T‘) 


rl  r2 


(2.15) 


(2.16) 


2  2  2  Q 

,  rl  r2  P  rl  r2  r2 

C3  *  “80  T~  0 

p  p 

Noting  that  Qr^,  the  axial  force  in  the  outer  cylinder,  can  be 


(2.16a) 


evaluated  from  equilibrium,  becomes  a  known  constant. 

The  problem  is  now  reduced  to  two  unknown  constants  and  which 

are  given  by  Equations  (2.14)  -  (2.15)  in  terms  of  the  redundant  moments  Mr 

and  Mr  . 

2p 


lp 


C.  Shell  Analysis 

Both  the  inner  and  outer  shells  have  an  identical  coordinate  system, 
pressure  loading  p,  rigidity  0C  and  modulus  E.  Hence,  the  governing  equation 
is  developed  for  a  typical  shell  with  radius  r. 


The  governing  4th  order  cylindrical  shell  equation,  giving  radial 


deflection  z  at  axial  position  y  is 

d4z  +  EV 

dy  r/D 


.  £_ 


(2. 


where 


Et 


0.  = 


'C  12(1  -  vZ) 

The  moment  and  shear  equations  for  a  cylindrical  shell  are 


v-°c4 

y  c  dy 


(2. 


d3z 


(2. 


dy 


P  rr 


The  particular  solution  to  Equation  (2.17)  is  z  *  -  ——  .  Let  $  * 


Et, 


4r  2D 
C 


be  a  shell  parameter.  Then,  Equation  (2.17)  reduces  to 


<2- 

dy  c 

The  solution  to  the  homogeneous  form  of  this  equation  is 

z  =  eBy(0^  cos  gy  +  Dg  sin  gy)  +  e'^Dg  cos  gy  +  sin  gy) 
where  D^t  D2,  D-j  and  D4  are  arbitrary  constants  of  integration. 

Replacing  the  exponential  functions  with  hyperbolic  functions  and 
adding  a  particular  solution,  the  solution  to  Equation  (2.17a)  becomes 
2 

prc 

z  -  -  —  +  D1  sin  (  gy )  sinh  (  gy )  +  02  sin  ( gy)  cosh  ( gy) 


+  03  cos  (gy)  sinh  (gy)  +  D4  cos  (gy)  cosh  (gy) 


(2 


By  symmetry,  if  y  is  changed  to  ( -y ) ,  the  solution  must  remain 
unchanged  i.e.  the  odd  functions  of  Equation  (2.20)  must  vanish.  Hence 


D, 


(2.21) 


2  3 

Subscript  1  refers  to  the  inner  shell  while  subscript  2  refers  to  the 
outer  shell.  The  four  constants,  two  from  the  inner  and  two  from  the  outer 


shell,  are  written  in  terms  of  the  redundant  moments  Mr  and  Mr  . 

lc  r2c 


Use  the  second  derivative  of  z  and  evaluate  moments  Mr  at  r  =  r  and 


lc 


1 


r  at  r  =  r  and  at  y  =  b. 


2c 


’rlc 


-  D„ 


d2z 


dy 


=  2Dc8iLAic°s(  Bib)cosh(  ^b)  -  B^sin(  B^b)sinh(  B^b)]  (2.22) 


V2c 


-  D 


d2z 


C  dy2 


2Dce2^l-A2c0S^  02b)cosh(  ^2b)  -  B2sin(  02b)sinh(  ^b)]  (2.23) 


D.  Compatibility 

Referring  to  the  free  body  sketch  of  Figure  (2.4),  we  can  write  the 
compatibility  conditions.  The  displacement  compatibility  is  expressed  as 

*1  =  0  at  y^  =  b  and  r  =  r^  (2.24) 

z2  ■  0  at  y2  =  b  and  r  =  r2  (2.25) 


The  slope  compatibility  is  given  by 

9cl  •  %>1  3t  r  ■  rl 
ec2  ■  ep2  at  r  =  r2 

The  moment  compatibility  is  given  by 


(2.26) 

(2.27) 


The  displacement  conditions  have  been  evaluated  using  Equation  (2.20) 
and  are  explicitly  shown  as  Equations  (2.30)  -  (2.31). 

prl2 

0  *  -  -ft —  +  A1sin($1b)  sinh  (^b)  +  Bj  cos  (Sjb)  cosh  (^b)  (2.30) 


0  *  -  ^ —  +  A2sin(  ^b)  sinh  (  ggb)  +  B2  cos  (B2b)  cosh  (  ^b)  (2.31) 

Explicit  results  for  slope  compatibility  are  as  follows.  Using  the 
first  derivative  of  Equations  (2.13)  and  (2.20)  evaluated  at  r  =  r  and  y  =  b, 
along  with  Equation  (2.26)  yields 

AlPlCsin(01b)cosh(  81b)  +  cos(  01b)sinh(  e^^b)]  +  B1e1[cos(  01b)sinh(  gjb) 

-  sin  ( 0jb)cosh(  0^b) ]  -  C2  (— )  -  (~)  =  -  +  ^3  (7—)  (2.32) 

1  1  p  1 

Using  the  first  derivative  of  Equations  (2.13)  and  (2.20)  evaluated  at 
r  *  r2  and  y  =  b  in  Equation  (2.27)  yields 

-A202[sin(82b)cosh(  +  cos(  B2b)sinh(  S2b ) ^“B2 02^-COS ^  e2b)sinh(  e2b^ 

^r2  1 

-  sin  (  &2b)cosh(  f^b) ]  -  C2( — 2 )  “  ^1  (7^) 


r?  p  r?  r? 

W  +  C3  2)[1  +  2ln  <77)] 

P  r,  1 


(2.33) 


Explicit  results  for  moment  compatibility  are  as  follows.  Equations 
(2.14)  and  (2.22)  evaluated  at  r  =  r^  using  Equation  (2.28)  yields 

2  2 
A1  C2Dc Bi  cos(01b)cosh(01b)]  -  B1[2DC01  si n (  01b)sinh(  Sjb)] 

20  0  r.2p  , 

-  C2(l  +  v)  — \  +  C^l  -  v)  -&2  =  -  (3  +  v)  +  C3  D  (i-^)  (2.34) 


.  •*  A.VV.%  \‘_V  1 


Equations  (2.15)  and  (2.23)  evaluated  at  r  ■  r2  using  Equation  (2.29) 


yields 

2  2 
A2[-2Dc02  cos(g2b)cosh(02b)]  +  B2[2Dce2  sin( 02b ) si nh (  ^b)] 

2D  D  r,2p 

+  Cg  (1  +  v)  — -  0^  (1  -  v)  -Ej  =  (3  +  v) 

rl  r2 


p 

"  CjDpC— -»r  ^  +  +  ~T  ^  +  v)^ 

r  i  r .  1 


(2.35) 


The  six  unknown  constants  of  integration  are  A^,  B^,  k^»  Bg,  and  C 
which  are  calculated  by  solving  Equations  (2.30)  -  (2.35)  simultaneously. 

In  summary,  given  the  system  parameters:  Fq,  p,  r^,  r 2>  tp,  tc>  E  an 
v,  we  calculate 


Qr2  =  (Fo  +  11  r22p) 

2  2  2  Q 

rl  r2  P  rl  r2  r2 

c3  - ?  T5 

p  p 


E  t. 


P  12(1-?) 

E  t  3 


Jr - T~ 

C  12(1- v) 


■  (-^M1/4 


*2  -  <-^>1/4 

4r2  Dc 


(2.36) 

(2.37) 

(2.38) 

(2.39) 

(2.40) 

(2.41) 


Solving  the  six  simultaneous  Equations  (2.30)  -  (2.35)  will  yield  the 
unknown  constants  of  integration.  The  deflection  of  the  plate  4  at  r  =  r,  is 


93 


then  computed  from  Equation  (2.13),  or 

4  4 


“W  ■  -  <r24  -  O  sfe- *  c3<f7>2 «"  «r> 
i  p  1  1 


+  C9  [(-i)2  -  1]  +  C,  tn  (-£) 


(2.42) 


For  a  unit  axial  load  P,  compute  the  deflection  a  using  Equation 
(2.42).  Then,  compute  the  effective  modulus  E'  of  the  equivalent  cylinder 
using  Equation  (2.2a). 


E.  Numerical  Results 


Let 


Define  a  nondimensional  bellows  parameter. 


X 


(2.43) 


1  x  Et 


(2.44) 


be  a  measure  of  system  response  i.e.  flexibility.  The  inverse  of  flexibility 
i.e.  stiffness  is 


(2.45) 


The  FORTRAN  program  SHELL. FOR,  whose  listing  is  given  in  Appendix  A, 
computes  the  effective  modulus  E'  for  a  given  set  of  system  parameters.  A 
numerical  example  is  shown  below.  Table  (2.1)  lists  the  numerical  values 
chosen  as  input  parameters. 

The  results  of  the  computation  are  as  follows: 


Design  curves  for  stiffness  S  vs.  nondimensional  parameter  X  are 
plotted  in  Figure  (2.6)  for  different  height  to  mean  radius  and  height  to 
one-fourth  length  of  corrugation  ratio.  These  curves  can  be  used  to  compute 
the  effective  modulus  E'  for  a  given  geometry  of  the  bellows. 


3.  ANALYSIS  OF  CANTILEVER  BEAM  ELEMENTS 


A.  Mathematical  Model 

A  typical  element  of  the  manipulator  arm  has  a  bellows  configuration 
with  the  ends  sealed  and  with  reinforcement  on  the  underside  as  shown  in 
Figure  (3.1).  To  achieve  the  bending  effect  of  the  manipulator  arm  in  either 
direction,  we  use  double  sided  elements  placed  back  to  back  as  shown  in  Figun 
(3.2). 

Each  element  of  the  arm  consists  of  several  corrugations,  and  the 
exact  number  is  determined  by  the  length  of  the  element  and  the  maximum 
desired  transverse  displacement  or  rotation  at  the  tip  of  the  cantilever. 

The  maximum  tip  rotation  will  typically  be  90°  or  180°  depending  on 
the  function  of  the  element.  Elements  that  are  designed  to  coil  around  the 
load  need  to  rotate  more  than  the  elements  that  help  in  the  lifting  of  the 
load. 

A  relationship  between  the  geometry,  length,  the  degree  of  rotation 
and  the  number  of  corrugations  in  the  double  sided  element  is  derived  as 
follows.  Assume  that  the  corrugations  on  the  convex  surface  just  touch  at  th 
maximum  degree  of  bend.  With  reference  to  Figure  (3.3),  these  two  relations 


Figure  3.2  Elevation  view  of  double-sided  bel 


forcement,  placed  back-to  back. 
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l  S 


/ 

v  .V 

i  'j 


L  *  R  ,4» 
r 


(R^  -  2r2)  ^  «  (2b  +  t)n 


(3.1) 

(3.2) 


where 


n  Is  the  number  of  corrugations 
b  is  the  one-fourth  length  of  the  corrugations 
t  is  the  maximum  degree  of  bend  of  the  element 
t  is  the  uniform  thickness  of  the  bellows 
L  is  the  length  of  the  element 
r^  is  the  outer  radius  of  the  bellows 
From  Equations  (3.1)  -  (3.2),  it  follows  that 

<1-  2r2>* 


(3.3) 


The  second  constraint  is  that  the  number  of  corrugations  and  the 


corrugation  length  should  satisfy 


(4b)n  =  L 


(3.4) 


The  two  unknowns  n,  b  are  solved  simultaneously  using  Equations  (3.3)- 


(3.4)  yielding 


7  -  2r2  * 


b 


(3.5) 


(3.6) 


The  elements  of  the  manipulator  arm  will  be  analyzed  as  a  cantilever 
beam  with  an  equivalent  modulus  E'  (refer  to  Chapter  II)  and  second  area 
moment  I  of  the  composite  single  or  double  sided  finger  element. 

The  composite  double  sided  element  is  symmetric  about  its  neutral 
axis.  Considering  half  of  this  section,  we  have  the  equivalent  cylinder,  slab 
and  reinforcement  as  shown  in  Figure  (3.1).  For  analysis  purpose,  we 
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consider  only  the  equival»~‘  cylinder,  a  flat  slab  of  modulus  E  and  two 
different  types  of  reinforcement:  n  threads  of  modulus  Er  and  radius  Rr,  or 
rectangular  reinforcement  of  modulus  Er,  as  shown  in  Figure  (3.4). 


For  double  sided  finger  elements  with  circular  threads,  the  second 
area  moment  is 

I  *  2{»R3te  +  2irRte(R  +  h  +  H)2  +  d^3  +  -  l)n(0.3927  Rp4)  }  (3.7a) 

For  double  sided  finger  elements  with  rectangular  reinforcements 

1  *  2{*R3te  +  2itRte(R  +  h  +  H)2  +  -J  d XH 3  +  -  1)  |  Sjhj3}  (3.7b) 

where  the  equivalent  thickness  of  the  corrugated  tube  is 


t 


e 


(3.8) 


B.  Deflection  Analysis 

The  manipulator  arm  is  made  up  of  a  number  of  elements  joined  end  to 
end.  The  function  of  the  arm  determines  the  size  and  number  of  elements  in 
the  arm.  The  cantilever  beam  analysis  of  a  typical  element  with  an  equivalent 
stiffness  El  is  presented. 

Large  deflection  is  taken  into  account.  Shear  deformations  are 
neglected.  The  shape  of  the  elastic  curve  is  found  from  the  exact 
differential  equation. 

This  cantilever  is  subjected  to  a  coplanar  terminal  loading  consisting 
of  an  axial  compressive  force  P$,  a  transverse  force  Pb  and  a  bending  moment 
Mb.  A  floating  coordinate  system  is  chosen  at  the  tip  of  the  cantilever  as 
shown  in  Figure  (3.5a).  In  previous  analysis  of  the  elastica  (Timoshenko  and 
Gere,  1961),  the  only  loading  was  an  axial  compressive  force  P*. 


/ 


MODULUS  t 
BASED  ON 
MODULUS  E 


SLAB  (MODULUS  E  ) 


The  exact  differential  equation  of  the  deflection  curve  is 


where 


(1  +  y'  )* 


EI  3? =  -  V  -  V  -  Mb 


Differentiating  Equation  (3.10)  with  respect  to  s  and  using  the 


relations 


we  obtain 


*  sin  0  =  cos  6 


2 

EI  |  *  -  P.  sin  6  -  PK  cos  9 
ds2  *  b 


Integrating  Equation  (3.11)  with  respect  to  ^  ds,  we  obtain 


7  E>  <5?>Z  *  %  9  -  pb  si"  9  +  Ki 


We  solve  for  the  constant  Kj  from  the  edge  condition:  At  s 
M. 

and  =  — -  .  Hence,  Equation  (3.12)  reduces  to 
as  EI 


M 

i  EI  (4r)2  =  P„(cos  0  -  cos  a)  +  P.  (sin  a  -  sin  0)  +  i 
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Since  ~  is  known  to  be  negative,  Equation  (3.13)  can  be  rearranged  as 

p  2  p  M 

ds  =  -  [2  — 4  (cos  0  -  cos  a)  +  — =1  (sin  a  -  sin  6)  +  (—■)^]"1^de  (3.13a) 
El  El  El 


Integrating,  nondimensionalizing  and  simplifying  Equation  (3.13a),  we 


obtain 


l  *  /  - 12 - , 

o  (A  cos  9  -  B  sin  0  +  c)' 


(3.14) 


where 


A  -  2  P* 

B  =2Pb 

C  ■  -2P i  cos  a  +  2  Pb  sin  a  +  Mb  2 


Let  A  cos  0  -  B  sin  0  =  D  cos  (0  +  0O).  Rearranging  Equation  (3.14), 


we  obtain 


l  -  /  - 12 - , 

o  [0  cos  {  0  +  0j  +  C] 


(3.15) 


Writing  0  +  Qq  =  2x  and  simplifying  Equation  (3.15)  we  obtain 


1 


a  +  0. 


(C  + 


T 


dx 


[1 


sin^x]1'2 


(3.16) 


where 

*2  _  2D_ 

*  "  C+D 

Equation  (3.16)  represents  an  indefinite  elliptic  integral  of  the 
first  kind.  Note  that  K  and  ^  depends  on  the  angle  a  and  the  loads. 

Solutions  of  this  elliptical  integral  are  available  in  tabulated  form,  given 
a,  for  pure  lateral,  longitudinal  or  moment  loading. 

Figure  (3.5b)  shows  an  element  subjected  to  a  pressure  loading  p  over 
a  uniform  area  A(=  ttR2)  acting  at  a  distance  d  and  a  dead  weight  W  acting  at 
an  angle  0. 

Resolving  these  loads  into  its  components,  we  obtain 

P*  =  -pA  cos  a  +  W  cos  0  (3.17a) 

Pb  =  pA  sin  a  -  W  sin  0  (3.17b) 

Mb  •  pAd  (3.17c) 

An  iterative  procedure  using  the  secant  method  and  the  IMSL  subroutine 
DCADRE  was  written  to  compute  the  unknown  angle  a  which  is  involved  in  the 
limits  of  Equation  (3.16).  A  FORTRAN  program  CANT,  whose  listing  is  given  in 
Appendix  B,  was  written  to  invoke  this  procedure  and  to  compute  the  shape  of 
the  elastic  curve.  The  initial  guess  value  of  o  is  taken  to  be  0. 

A  numerical  example  for  a  typical  element  was  computed.  The  length 
and  stiffness  of  the  element  were  chosen  as  l  =  1  in  and  El  =  1  lb  in^, 
respectively.  The  pressure  loading  acts  at  an  eccentricity  of  0.1  times  the 
length.  For  increasing  pressure  loading,  the  corresponding  values  of  the 
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eccentric  end  load  P  and  the  end  moment  chosen  are  listed  in  Table  (3.1). 
The  calculated  deflection  curves  at  these  loads  are  shown  in  Figure  (3.6). 

The  deflection  pattern  resembles  that  of  the  elastica.  Table  (3.1)  shows  the 
tip  angle  o  corresponding  to  each  level  of  loading. 


Eccentric 
end  load 

P  lb 

End 

moment 

Mb  lb-in 

Calculated 
tip  angle 
a  deg. 

0 

0 

0 

2 

0.2 

16 

4 

0.4 

41 

6 

0.6 

78 

8 

0.8 

124 

10 

1.0 

173 

Table  3.1  End  loads  and  the  corresponding  tip  angle  for 
W  *  0,  i  *  1  in,  El  *  1  lb-in2  and  d /*,  =  0.1. 


4.  ELEMENT  STRINGS 


A.  Transformation  Equations 

The  manipulator  arm  is  made  up  of  a  number  of  elements  connected  by 
rigid  links.  One  end  of  the  arm  is  the  supporting  base  of  the  cantilever  beam 
while  the  last  few  elements  near  the  tip  act  as  a  gripper  or  end  effector. 
These  tip  elements  grip  the  payload  by  coiling  around  it.  The  points  of 
contact  between  the  element  and  the  load  are  sufficient  to  prevent  the  load 
from  slipping.  The  rigid  links  do  not  bend.  They  merely  translate  and 
rotate. 

Let  the  number  of  elements  in  the  "lifting"  arm  be  n  and  the  elements 
numbered  consecutively  from  the  supporting  base  l,2...n.  Each  element  of  the 
arm  is  modeled  as  a  cantilever  beam  starting  from  the  base.  The  end  of  the 
rigid  link  after  element  i  is  the  fixed  end  of  the  cantilever  beam  for  element 
i  +  1.  Given  a  coplanar  terminal  loading,  we  can  calculate  the  end 
deflections  of  that  element  by  the  methods  presented  in  the  last  chapter. 
Hence,  we  have  determined  the  end  coordinates  of  the  element  with  respect  to 
the  local  coordinate  system. 

Determining  the  position  and  orientation  of  any  element  in  the  arm 
relative  to  its  base  coordinate  requires  the  transformation  of  coordinates 


through  all  other  elements  and  links  between  the  base  reference  and  the 
element  whose  coordinates  are  being  determined.  Similarly,  by  carrying  out 
the  transformation  for  each  of  the  elements,  we  get  the  overall  orientation  of 
the  "lifting"  arm  in  the  global  coordinate  system.  Figure  (4.1)  shows  a 
relation  between  the  local  and  global  coordinate  system. 

A  transformation  matrix  relating  the  (i)  frame  to  the  (i+1)  frame  is 
of  the  form: 

X,tll  rcos(vi>  s(n<Vi>  xm,p 

+ 

Yi+1  _sin(o^+j)  cos(aj+i)  Yi+ifp 

where 

ai+^  is  the  angle  measured  clockwise  from  frame  (i)  to 

frame  (i+1). 

x-j+l»Yi+l  are  the  coordinates  of  the  origin  of  frame  (i+1)  in  the 
coordinate  system  of  frame  (i). 

*i+l,p»Yi+l,p  are  the  coordinates  of  any  point  P  in  frame  (i+1). 

Xi  n»Yi  n  are  coordinates  of  the  same  point  P  in  frame  (i). 

A  FORTRAN  program  TRANS,  whose  listing  is  given  in  Appendix  C,  models 
each  lifting  element  as  a  cantilever  beam  in  the  local  coordinate  system.  The 
slope  computed  at  the  tip  of  element  i  is  the  angle  a.+^  which  is  used  to 
transform  the  element  (i+1)  from  its  local  coordinate  system  to  that  of 
element  i.  The  position  of  the  tip  of  the  element  in  the  local  coordinate 
system  is  then  computed.  Provision  is  made  for  the  transformation  between  an 
element  and  the  link  before  it.  As  the  link  is  rigid,  the  angle  a.  ,  is  zero 
degrees.  The  transformation  matrix  relating  the  tip  of  element  i  and  the  base 
reference  is  the  product  of  all  matrices  relating  the  elements  and  the  links 


FIRST  ELEMENT 


NEUTRAL  AXIS 


(a) 


RIGID 

LINK 


Figure  4.1  First  element  showing  coincidence  of  local  and  global  axes 


(a);  rigid  link  showing  translation  and  second  element  show¬ 
ing  displaced  local  axis  (b). 


between  them.  Hence,  the  program  TRANS  determines  the  position  of  the  tip  of 
each  lifting  element  of  the  arm  in  the  global  coordinate  system. 

B.  Typical  Design 

Figure  (4.2)  shows  the  motion  of  a  manipulator  arm.  Double  or  single 
sided  elements  could  be  used  as  the  structural  section  of  the  elements.  The 
stiffness  of  double  sided  elements  were  derived  in  part  A  of  Chapter  III.  The 
advantage  of  double  sided  elements  is  that  they  may  bend  the  arm  in  either  the 
clockwise  or  counter  clockwise  direction.  For  instance,  if  the  elements  on 
the  left  side  are  pressurized,  the  arm  will  bend  in  the  counter  clockwise 
di recti  on. 

Given  a  destination  coordinate  in  the  feasible  work  space,  we  want  to 
know  a  pressure  history  in  all  the  elements  that  will  move  the  load  from  the 
starting  to  the  given  destination  coordinate.  Note  that  there  are  many 
possible  solutions  or  different  pressure  histories  that  could  achieve  the  same 
result.  The  solution  procedure  described  below  picks  a  solution  that  may  not 
necessarily  be  the  one  where  the  payload  tranverses  the  least  path  distance  to 
achieve  its  target  point.  Mahajan  (1985)  investigated  the  problem  of  optimal 
paths. 

The  end  effector  is  designed  separately.  The  number  of  these  elements 
is  determined  by  the  length  of  the  element  and  the  size  of  the  object  to  be 
lifted.  The  end  effector  bends  in  the  opposite  direction  to  that  of  the 
lifting  arm.  The  angle  0i  to  which  these  end  effector  elements  should  coil 
around  the  load  to  prevent  it  from  slipping  is  determined  by  a  simple 
geometrical  drawing.  The  slope  a  at  the  end  of  each  of  these  elements  is 


LOWER 

ELEMENTS 


where  Ne  is  the  number  of  end  effector  elements.  Then,  using  the  program 
CANT,  we  make  runs  of  varying  pressure  levels.  Since  we  do  not  know  the  point 
of  contact  between  the  load  and  the  end  effector  elements,  we  assume  the  worst 
case  or  the  case  of  highest  moment.  That  is 

0  =*  180°  -  a  (4.2) 

M0  =  W(R*  +  2R  +  d  +  H)  (4.3a) 

M0  =  W  R£  (4.3b) 

The  moment  expressions  in  Equations  (4.3a)  -  (4.3b)  are  for  the  double  and 
single  sided  elements  respectively.  Hence,  the  pressure  history  in  the  end 
effector  elements  is  determined  by  that  pressure  level  whose  slope  at  the  tip 
of  the  element  is  closest  to  a. 

The  lifting  elements  are  designed  as  an  element  string.  Initially, 
the  moment  due  to  the  payload  is  assumed  to  be  zero.  Referring  to  Figure 
(3.5b),  the  angle  0  at  which  the  load  acts  is  updated  from  the  slopes  of  all 
the  previous  elements.  For  the  ith  element,  the  angle  0  is  given  by 


i-1 

0i  »  180°  -  I  « 
k=l  * 

i*2,3...n 

(4.4a) 

0j  =  180° 

(4.4b) 

Then,  for  a  given  set  of  internal  pressures  in  each  of  the  elements,  we 
compute  the  position  of  the  tip  of  the  element  in  the  global  coordinate  system 
using  the  transformation  matrices  developed  in  part  A.  Knowing  the  geometry 
of  the  link  joining  the  lifting  elements  to  the  end  effector  and  the  radius  of 
the  load,  we  can  compute  the  coordinates  of  the  load.  Referring  to  Figure 
(4.3),  we  compute  the  coordinates  of  the  load  as  follows. 

RTPL  =  2  Rb  +  H  +  Rl  (4.5) 

-1  XL 

*  =  tan 


lV.V. 


.  > 


(4.6) 
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n 

Gx  =  x-coordinate  of  nth  element  +  Rjpl  sin(  £  a  +  $)  (4.7) 

k-1  k 

n 

Gy  *  y-coordinate  of  nth  element  -  Rjpl  cos(  l  a  +  $)  (4.8) 

k-1  k 

Knowing  the  coordinates  of  the  load,  we  can  compute  the  moment  due  to  the  load 
acting  on  each  element.  The  entire  procedure  is  repeated.  The  following 
inequalities  must  be  simultaneously  satisfied  for  the  specified  tolerance 
limit  5. 

|GxCurrent  -  G*earlierl  <  5  (4.9) 

|Gycurrent  "  Gyearlierl  K  5  (4.10) 

The  iteration  stops  when  both  the  inequalities  are  satisfied.  Hence,  for  a 
given  set  of  pressures,  we  can  compute  the  global  element  deflections  of  the 
manipulator  arm. 

A  systematic  procedure  is  adopted  to  vary  the  pressure  in  the  elements 
to  move  the  load  from  the  initial  or  present  coordinates  of  the  load  (Xp,  Yp) 
to  the  destination  coordinates  (Xf,  Yf).  The  present  position  of  the  load  can 
be  in  any  one  of  the  four  regions,  as  shown  in  Figure  (4.4).  With  reference 
to  Figure  (4.2). 

i)  The  load  is  in  the  upper-right  region,  or 
Xp  <  Xf 
Yp  >  Yf 

Then  decrease  the  pressure  in  all  the  elements, 
ii)  The  load  is  in  the  upper-left  region,  or 
Xp  <  Xf 
Yp  <  Yf 

Then  decrease  the  pressure  in  the  lower  elements. 


1 


0 


UPPER -LEFT 
REGION 


LOWER-LEFT 

REGION 


Y 


i 

UPPER  -  RIGHT 
REGION 

(Xf,Yf)  FINAL  LOAD  POSITION 

LOWER-RIGHT 
REGION 


(Xp.Yp)  • 

PRESENT  LOAD  POSITION 


Figure  4.4  The  four  possible  regions  in  which  the  load  could  be 
presently  positioned. 


iii)  The  load  is  in  the  lower-right  region,  or 
Xp  >  Xf 
Yp  >  Yf 

Then  increase  the  pressure  in  the  lower  elements, 
iv)  The  load  is  in  the  lower-left  region,  or 
Xp  >  Xf 
YP  <  Yf 

Then  increase  the  pressure  in  all  the  elements. 

The  iteration  procedure  is  repeated  until  the  following  conditions  are 
satisfied: 

|Xp  -  Xf |  <  Cl  (4.11) 

I Yp  -  Yf |  <  c2  (4.12) 

where  ei  and  E2  are  specified  tolerance  limits. 

C.  Numerical  Results 

The  manipulator  arm  shown  in  Figure  (4.2)  is  composed  of  six  single¬ 
sided  SIMRIT  finger  elements  connected  by  rigid  links.  The  first  four 
elements  are  the  lifting  elements  while  the  last  two  form  the  end  effector. 
There  are  three  different  types  of  elements  A,  B  and  C.  Figure  (4.5)  shows  a 
cross  section  through  the  SIMRIT  finger.  Table  (4.1)  summarizes  the  geometry 
for  these  elements.  The  object  to  be  lifted  is  a  cylinder  of  radius  1.25  in. 
and  weighing  0.5  lb.  The  stiffness  EI  for  these  elements  were  derived  by 
Wilson  (1986)  and  calculated  by  Walker  (1986). 

From  a  scale  drawing,  refer  to  Figure  (4.3),  we  note  that  the  tip 
angle  relative  to  its  own  base  for  the  last  2  elements  (type  C)  is  206°.  In 
addition  to  the  geometry  of  this  element,  the  other  input  parameters  to  the 
program  CANT  are 


Figure  4.5  Cross  section  through  a  SIMRIT  finger  showing  the 


geometry 


Dimension  Size  A  Size  B  Size  C 

(unit) 


b  (in) 

0.1150 

0.0828 

0.0603 

dl  (in) 

1.563 

1.063 

0.8000 

d2  (in) 

0.7087 

0.4930 

0.3260 

d3  (in) 

0.2314 

0.1882 

0.1380 

d4  (in) 

0.0948 

0.0694 

0.0700 

h  (in) 

0.1745 

0.1286 

0.0960 

L  (in) 

3.875 

2.688 

1.730 

R  (in) 

0.6875 

0.4625 

0.2550 

t  (in) 

0.0943 

0.0826 

0.0700 

a  (rad) 

2.950 

2.827 

3.013 

yp  (in) 

0.7156 

0.4875 

0.3133 

Ap  (in2) 

1.296 

0.5366 

0.2150 

I  (in4) 

0.3700 

0.0730 

0.0150 

E'  (lb/in2) 

137.2 

189.7 

240.9 

E  (lb/in2) 

2120. 

2120. 

2120. 

El  (lb-in2) 

50.77 

13.92 

3.640 

Table  4.1  Geometry  for  three  different  types  of  SIMRIT  finger  elements 
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3  =  77° 

M0  =  -0.656  lb-in 
W  =  0.5  lb. 

The  pressure  loading  that  corresponds  to  a  maximum  slope  of  103°  is  37.2  psi. 

The  lifting  elements  are  analyzed  using  the  program  TRANS.  The  first 
two  elements  are  of  type  A  while  the  other  two  are  of  type  B.  Initially,  the 
moment  due  to  the  load  is  assumed  to  be  zero  in  all  of  the  elements.  Noting 
that  the  rigid  link  does  not  bend  (i.e.  a=  0°),  the  transformation  matrix  is 
given  as  input  data.  Each  of  the  other  elements  is  modeled  as  a  cantilever 
beam  in  its  local  axis  and  the  corresponding  transformation  matrix  is 
computed.  The  subroutine  TRANS  is  then  invoked  to  obtain  the  coordinates  of 
the  tip  of  each  lifting  element  in  the  global  coordinate  system.  Then,  using 
Equations  (4.5)  -  (4.8),  the  coordinates  of  the  center  of  gravity  of  the  load 
is  computed.  The  moment  due  to  the  load  is  then  updated.  The  iteration 
continues  till  the  conditions  given  in  Equations  (4.9)  -  (4.10)  are  satisfied. 
In  addition  to  the  geometry  of  the  two  elements,  the  other  input  parameters  to 
the  program  TRANS  are 

Rl  =  1.25  in. 

X|_  =  0.5  in. 

YL  =  0.75  in. 

H  =  0.1  in. 

The  results  of  the  computation  is  summarized  in  Table  (4.2).  The  trajectories 
corresponding  to  the  pressure  loading  shown  in  Table  (4.2)  are  plotted  in 
Figure  (4.6). 

The  program  TRANS  can  handle  many  different  options.  An  element  in 
the  arm  may  or  may  not  be  pressurized.  The  elements  can  be  of  varying 
geometry.  However,  for  large  nondimensional  pressure  parameters  (Pj,  Pb  >  10) 


cood. of  element 


when  subjected  to  different  pressures. 
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the  results  of  this  program  should  be  used  cautiously  as  numerical  instability 
may  occur.  This  needs  further  investigation. 

Using  the  systematic  procedure  described  in  part  B,  one  possible 
pressure  loading  to  lift  the  load  to  a  vertical  position  (-0.3,  7.2)  is 


Element 


Pressure  (psi) 


37.20 

37.20 


•JjlJ  *  jl.  'rf-  - 
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b.  CONCLUSIONS  AND  DISCUSSION 


The  design  curves  for  stiffness  vs_.  the  nondimensional  parameter  A  was 
found  to  be  within  ±  6%  of  that  predicted  by  the  Haringx  model.  In  his  paper, 
Wilson  (1984-b),  has  shown  that  the  Haringx  model  correlates  well  with  that  of 
the  experimental  data. 

The  elastic  curve  of  the  deflection  pattern  of  the  cantilever 
subjected  to  an  increasing  pressure  loading  was  found  to  resemble  that  of  the 
elastica.  The  elastica  problem  (Timoshenko  and  Gere,  1961),  in  which  only  the 
longitudinal  loads  were  considered,  were  verified  using  this  model. 

The  calculated  trajectories  of  a  six  SIMRIT  finger  arm  subjected  to 
increasing  pressure  histories  are  shown.  These  trajectories  represent 
feasible  equilibrium  positions  of  the  arm.  The  actual  position  of  the  arm  in 
the  global  coordinate  system  will  be  verified  by  experiment  in  due  course. 

There  are  several  possible  extensions  of  this  research.  The  numerical 
example  shown  in  Chapter  IV  consists  of  single  sided  elements  with  the  neutral 
axis  for  all  the  lifting  elements  on  the  same  side.  In  this  configuration, 
the  arm  can  move  only  in  the  counter  clockwise  direction  through  180°  in  the 
same  plane.  To  achieve  reversal  of  curvature,  we  could  either  pressurize 
double  sided  elements  judiciously  or  have  the  neutral  axis  for  part  of  the  arm 
on  one  side  and  for  the  remaining  part  of  the  arm  on  the  opposite  side.  With 
the  reversal  of  curvature,  the  arm  could  move  ±  180°  in  the  same  plane.  The 
mathematical  tTPdel  for  plane  motion  can  be  extended  to  include  out-of-plane 
motion.  This  can  be  achieved  by  arranging  the  bellows  in  a  satellite 
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configuration  as  shown  in  Figure  (5.1).  In  the  case  of  an  even  number  of  the 
bellows  (Figure  (5.1a)),  we  can  pressurize  the  bellows  (1,2, 3, 4)  to  achieve 
bending  about  the  axis  A-A.  Hence,  the  flexible  arm  can  have  a  work  space 
consisting  of  volume  of  revolutions  in  each  of  the  four  quadrants. 

Provision  is  made  in  the  program  TRANS  to  include  the  self  weight  of 
the  element.  The  self  weight  of  the  element  is  assumed  to  be  a  concentrated 
mass  acting  in  the  direction  of  gravity  at  the  tip  of  the  element. 

A  'Learning  Program*  was  developed  by  (Mahajan,  1985)  to  vary  the 
pressure  in  a  systematic  way  to  achieve  rapid  convergence  to  the  destination 
coordinate.  In  his  model,  only  the  end  moment  was  considered.  Mahajan' s  pro¬ 
gram  could  be  extended  to  include  the  mathematical  model  presented  in  Chapter 
IV,  a  model  that  includes  end  moment,  longitudinal  and  transverse  loading. 

Computational  efficiency  could  be  improved.  The  numerical  instability 
that  occurs  in  the  program  TRANS  for  nondimensional  pressure  loading  of'P *  > 

10  or  Pfc  >  10  could  be  overcome.  In  addition,  the  number  of  iterations 
required  to  compute  the  tip  angle  a  in  Equation  (3.16)  may  be  reduced.  We 
could  form  a  data  base  for  the  three  different  nondimensional  loading 
parameters  P^,  P^  and  Mj).  The  angle  a  should  lie  between  0°  and  180°.  For  a 
particular  combination  of  Pz,  Pt,  and  Mb,  we  can  store  that  value  of  a  which 
satisfies  Equation  (3.16).  Hence,  the  program  TRANS  can  pick  the  right  angle 
a  from  this  data  base  without  having  to  perform  the  numerous  iterations. 

Another  bending  geometry  that  could  be  investigated  is  that  of  a 
tapered  section,  with  the  broader  end  at  the  supporting  base.  This  could  be 
an  optimal  bending  geometry  since  the  moment  arm  due  to  the  load  will  be 
higher  at  the  supporting  base,  decreasing  towards  the  gripper.  Such  tapered 
configurations  would  have  a  much  higher  payload  to  self  weight  ratio  than 
uniform  configurations. 
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APPENDIX  A 

COMPUTER  PROGRAM  FOR  THE  COMPUTATION  OF  REDUCED  MODULUS 


r»ea  1  atlC,  10).  *<10>,y  UO),  mu, -irie.mrEe, —IP,  fr£b 

c  <  10,  10)  .  f .  eo*.  cm**,  temp  (10)  .  w.  iamoc* 
int«p»i'  *(10> 

Inout  and  out  out  file  o»ci»r»tion*. 

OPEN  <7, PILE-’ SHE.  DOT’  ) 

OPEN (8, FILE-’ SHE.  OUT’  ) 

Read  in  axial  load  and  internal  pressure  on  pel  low*. 

read  <7, 1 ) FO,  d 
writ*<8, 3) c0,  p 

format (//*  Pxial  load-’ , F10. 4/’  Pressure-’ , FlO. 4//) 

Raad  in  Young’*  Modulu*  and  Poiwen’*  ratio, 
read (7, 1 ) E, mu 

writ#<8, £1)E, mu  _  A  .. 

format (’  Young*  ftodul u.-’ FlO. 4/’  Poi»*. rat  10-  ,F10.4> 

Raad  in  numbar  of  interval*  along  Plata,  inner  and  outar 
cylinder*  at  which  deflection*  are  cesirec. 

read  (7,  22)  int 

wr  it  *  ( 8,  t3)  int  .  .  ta//) 

format (/*  Interval  at  which  oef.n.  l*  de*irad  , I 

Read  in  mean  radius. 

raad  <7,  1 )  Rm 
write (8,  43) Rm 

formate/*  Mean  rad iu*-’  , FlO.  4//) 

Loop  to  control  variation,  in  depth  of  corrugation, 
do  100  kk-1,2 

Read  in  half  deoth  of  corrugation. 

read (7, 1 )  h 
r 1— Rm— h 
r&— Rmrh 

writ*  (8, 4)  rl,  r2  ..  ,  _iri  .. 

format (/*  Inner  cyl.  rad i u*-’  ,  FlO.  4/’  Outer  cyl.  radius-  ,F10.4) 

Looo  to  control  variation*  of  corrugation  length. 

do  101  JJ-1,5 

Read  in  one-forth  length  of  corrugat ion. 

read  <7, 1 ) b 
write (6, 10) o 

format (/’  Half  lenth  of  cy 1 .  -’  . F10. 4) 

Thickness  of  plat*  and  cylmoers  is  as*um*o  to  be  the  same. 
loop  to  control  variations  in  cylinder  thickness. 
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doer  rtnt 


do  i02  1 i *  1 . 7 

freed  in  cylinder  cmcknees.  -  P|  -v*!-dble  to  DTIC  doe*  nee 

*“**  Ml  build.  ,.Koduc,  J"* 

rtaa i 7 , 1 > tc 
to*te 

writeC8.20)tp.tc 

formati/'  Thick.  oi  Diaie=  '  .  FIO  .  •*/ '  Thick,  ot  she  1 1  ■  '  .  F10  .  4  ) 

format.  ( 2014 ) 

format(bF10.4> 


Number  of  unknown  constants  of  integration  is  6. 
n*6 

rstlo*r2/r 1 
do  12  i*l , n 
do  12  3*1. n 
A  < 1 , 3  > "0 . 0 

Compute  the  known  constant  and  other  parameters  to  be  used 
in  the  formulation  of  the  linear  system. 

0r2* (FO-3. 14159*r2*r2*p> / <2*3. I4lb9*r2> 
Dp=E*tp*tp»tp/(l2»<l-mu*mu> ) 

Dc*E* tc* tc* tc/ <12* (l-mu*mu)  ) 
e3irl«rl*r2/Dp*<  r2*p/8-Qr 2/4 > 

Bl*sqrt  isqrUE»tc/  (4»rl*rl*Dc>  >  ) 
B2“sqrt<sqrt<E«tc/<4»r2*r2*Dc> ) ) 

Formulation  of  the  6x6  simultaneous  linear  system. 

y(l)*p*rl«rl/(E*tc) 

y(2)*p*r2»r2/(E*tc> 

templ*-rl*rl*rl»p/16 

y(3>*templ/Dp*c3/rl 

y<4)*(3*mu>»(templ/rl*Dp*c3/(rl*rl)> 

temp2*-r2*r2*r2*p/16 

y (5) *temp2/Dp*c3*r2/ <  r 1 »r 1 ) *<l*2*log(r2/rl>> 
temp3  =  r2*r2*p*  <3*mu>/16 

temp4=c3*Dp/(rl*rl>*((3*mu)*2*il*mu)* log  <  r2/r 1 ) ) 

y<6>*temp3-temp4 

bbl=Bl *b 

bb2=B2*b 

a<l,l>*ain(bbl>*slnh(bbl> 
a(l,2)*cos<tbl)*coah<bbl) 
a(2,3>*sin(bb2>»slnh(bb2> 
a(2,4)=coa(bb2> *cosh <  bb2 ) 

a<3,l>=bl*<ain(bbl)»coah(bbl>»coa(bbl)»sinh<bbl>> 

a (3, 2) *B1* (cos(bbl) *slnh  <  bbl ) -sin  <bbl ) *cosh(bbl ) ) 

a(3,5>*-l/rl 

a(3,6>*-2/rl 

a<4,l)*2*Dc*Bl*bl »cos ( bbl ) -cosh  C bbl ) 
a  <4 , 2) *-2*Dc*Bl »Bl*ain (bbl )»sinn(bbl) 
a(4,5)*0p/(rl*rl)*(l-mu) 
a(4,6)*-2*Dp/(rl*rl)*(l*mu> 

a(5,3)*-B2«<sinlbb2) »cosh (bb2) -cos (bb2) *sinh(bb2) > 
a (5 . 4 ) * -B2*  <cos( bb2; - sinn ( bb2 ) -  sin  <  bb2 ) *cosh ( bb2 ) ) 
a (5 . 5 ) *  - 1 /r 2 
a(5,6>=-2*r2/(rl»rii 

a(6.3)=-2*Dc*B2»B2*cos(bb2) *cosh l bb2 ) 
a<6.4)=2*Dc*32*B2*sin<bb2)»sinh(bb2) 
a(b,5)*-Dp/(r2*r2)*il-mu> 
a(6.6)*2*Dp/(rl*rli*(l-muJ 


do  40  i *  1 . n 
s  ( 1  ;  -  i 


Copy  available  to  t 

pc.:  ii  i  1 


/  ■ 


1  O'? 


E 


x  v  i  > *y i i ) 
co  40  t * i . n 

40  c!l,j>«e!i,j) 

eps* le- 15 

c  Call  subroutine*  lr  and  ib  to  solve  the  b  x  b  simultaneous  ecuations. 

call  1 r < c , n , a , eps , aucess )  • 

call  j'Mc,n,s,x) 

c  Re-order  solution  vector  x. 

do  47  l *  1 , n 

47  temp < i ) *x t a ( i ) ) 
do  48  l»l,n 

48  X  !  1 ) ■ temp ! 1 ) 

wn te  !8, 6) lx  1 1  > , 1» 1 ,n) 

6  for*at(//'  The  constants  A1 . B1 , A2 . B2 , Cl , C2  ere  resp . : ' //6cl3 . 6) 

c  Compute  slope  and  moment  at  junctions  between  plate  and 

c  inner  and  outer  cylinders. 

Thetscl«x(l)*a(3,l)*x(2)»a(3,2) 

Thetac2«x (3) *o (5.3) *x  14)  *a  <5 , 4  > 

Thetapl *2/rl *x  6 >  *x ! 5 > /r 1 - y ! 3 ) 

Thetep2«2»r2*x  !b)/!rl«rl)*x!5)/r2*y!5) 

Kr2c»x!3)*a!6,3)*x!4)«e!6,4) 

Mrlc«x!l)*a!4,l)*x!2)*al4,2) 

Mrlp=y!4)-a!4,6)»x!6)-a!4,5)*x!5) 

Mr2p«y !6) -a !6,6) *x (b) -a  t6, 5) «x (5) 
write!8,25)Krlc 

25  format!//,'  Moment  Mrlc  at  r  =  rl  of  inner  cy  1  inder '  ,  FI  0 . 4  ) 
wrlte!8,26)Mr2c 

26  formate/,'  Moment  Mr2c  at  r«r2  of  outer  cy 1 inder ' , F10 . 4 ) 
write(8,27)Thetacl 

27  formate/,'  slope  Thetacl  at  r=rl  of  inner  cy 1 inder ' , F10 . 4 ) 
write<8,28)Thetac2 

28  format!/,'  slope  Thetac2  at  r«r2  of  outer  cy 1 inder ' , F10 . 4 ) 
write <8, 29) Thetapl 

29  format!/,'  slope  Thetapl  at  r*rl  of  plate ' , F10 . 4 > 
write!8,30)Thetap2 

30  format!/,'  slope  Thetap2  at  r«r2  of  plate ' , F10 . 4 ) 
write!8,31)Mrlp 

31  format!/,'  Moment  Mrlp  at  r-rl  of  plate ' , F10 . 4 ) 
write !8, 32)Mr2p 

32  format!/,'  Moment  Kr2p  at  r*r2  of  plate ' , FI 0 . 4 ) 
write !8, 35) 

35  format!//'  radii  deflection'//) 

hh« (r2-rl ) /int 
r*r  1 

c  Compute  axial  deflections  of  the  plate  at  the  specified  intervals. 

do  33  i = 1 , int* 1 
wl«-!r»*4-rl«*4)*p/!b4*Dp) 
rr=r/r 1 

w2*c3*rr*rr » log !rr) 
w3= ! rr*rr- 1 ) *x!6) 
w4=ioa!rr)*x!5) 
w= wl - w2 * w3 * w4 
wr i te ! 8 , 34 ) r , w 

33  r=r-hh 
del ta  »w 
write ! 8 , 37 ) 

37  format!//'  Outer  cyj  dist  deflection''//) 

Copy  available  to  DTJC  does  not 
.  permit  fully  legible  reproduction 
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nn— p/ int 

C  ist-C. 

cc.nst«-s*r2#f£/  (;+tC) 

z  Compute  oef lections  of  outer  cy.inoe-  at  fit  istcfitc  lrte-va.s. 

ao  36 

032-62*01 st 

2-eonst*x (3) *sin  <  0  02) •sinft <  002  >  *x (4  >  *cos ( dd2) *eosn  <  dd2 ) 
write (8,  34) diet,  z 
36  dist-dist*hh 

dist-0. 
write  <8, 38) 

38  format (//’  Inner  cyl  dist  deflection’//) 
const— p*rl*rl/  <£*tc> 

c  Compute  deflections  of  inner  cylinder  at  tne  specifiec  intervals. 

do  39  i«l, int*l 
bbl-Bl*dist 

r*con*t*x < 1 ) *sin  ( obi ) *si nh  <bbl >  *x (2) *cos ( bbl ) *eosn ( doi ) 
write  <8,  34) dist,  z 

39  dist— dist*hh 

c  Compute  effective  modulus  and  modulus  ratio. 

El—  <F0*3.  14157*Rm*Rm*p) *b/  (3.  14157*Rm*to*oelta> 

EEl— E/El 

write  <8, 41)E1, EEl 

41  format <//’  Effective  Modulus  El-’ ,  dl3.  6//’  Modulus  ratio  EEl**, 

1  dl3. 6//) 

34  format (2d 13. 6/) 

c  Compute  nondimensional  parameter  lambda,  flexibility  and  stiffness. 

lambda— tc/ (b*b) *Rm 

f 1 ex-EE l*te/Rm 

stiff— 1 /flex 

write <8,  49) stiff, lambda 

49  formate//,’  STIFFNESS-’ , F15. 6/’  LOMBDft— ’ , F15. 6//) 

102  continue 

101  continue 

100  continue 

stop 
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THIS  SUBROUTINE  RAC”ORISSS  THE  COEFFICIENT  - 

IN-C  *~S  .-R  e«C*ORS  AND  WRITES  *hE»  OVER  MATRIX  A. 
3RER  :nd:=A-ES  OPTIONS:  SREC-SUJ-C  rOR  Nr  AOR-'IA-  = 
PREC-1  FQR  PARTIAL  PIVOT  SEARCH. 


:  vc 


tuDi*cutint  in  (  a.  n,  a,  *oa,  xuc#»* ) 
intagi’*  »(10)  (  t««iPi  praf 
Character  luctii*l 
eiaanaion  aCO,  10) 


or**f»l 

if («<1>. ao.O)than 
or***1 (  1  > 
•<1>-1 
•no  if 


THIS  BLOCK  INTERCHANGES  ROWS  IN  THE  CASE  OF  PARTIAL  PI 


VO"  Ss.SSCh. 


10 


if  (oraf.  na.  0)  than 

IMX*0. 

do  10  i»l,  n 

if  <«M  (a  <  i,  1)  >.  g*.  «m»K )  than 
«M»>abi  (a(i|  1  > ) 
k»i 
and  if 
cent inui 
if  (k.na.  llthan 
(  1  ) 

*(l)-«(k) 
s (k)-tamD 
•nd  i  f 
and  if 


PIVOTING  ABOUT  PIVOT  ELEMENT  (1,1) 


do  30  i*£, n 
n£~* ( i ) 

if (aba <a(»<l >, 1) ).  It.  ap*) than 

aueaaa*'  n’ 

go  to  130 


1 V 

a(n£, l>-a(n£, l)/a<S<l>, 1> 

•ndif 

Si 

M 

30 

rontinua 

v 

c 

PIVOTING  ABOUT  ANY  GENERAL  PIVOT  ELEMENT 

% 

do  40  kl“£. n 

do  50  k£-kl,n 

*um»0. 0 

do  60  j»l, kl-1 

60 

»ur»**ui»*a  <»  <k£>  ,  J  >  *a  (*  <  J  >  ,  Kl  > 

3S 

50 

a<s(k£),  kl)-a(*<k£) ,  ki)-*um 
if (Ml. It. n) tnan 

if  (Dnaf .  na.  0>th*n 

aaiaxM). 

. • 

t.v 

do  70  n3*kl,n 

if (abs (a (* (k3) , kl ) ) .  pa. amax) than 
aaaxMb*  (a  (s  ( k2 ) ,  kl )  ) 

(<4»h3 
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painil  hilly  1 «0»W»  todac»tm 
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30 


:oo 

30 


*0 


•nc  \f 
eont 1 nu» 

if (k4.  n«.  kl>tn»n 

tt«D*l (HI  ) 

•  (  «1  >  ••  (•«*) 
s  (n«)ntMp 
•rid  if 
•rid  if 

if  <SDS<*<»<kl>,  kl)  >.  it.»0*)tn»r, 

fUC*fl>'  «’ 

po  to  130 

•1** 

oo  80  il-kl*l,n 

*<s(il),kl>-*(s<il>»kl)/*l*<kl>’kl) 
do  90  jl»kl*l,n 

sum— 0. 

do  100  i£-l,kl-l 

sum«sum*«(s<kl>, i2> *• <s ( i£) ,  J 1 > 

* (S  <kl > , J 1 > — • (s(kl),jl>  —sum 
•no  if 
•ndif 
eont  i  nut 


:%■ 
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-MIS  SUBROUTINE  DOES  *ORwPRD  aro  BSCKWORD  OOSSES  On 

ONY  B  VECTOR.  THE  RESJ--S  0-  ”!'HE  DOSSES  ORE  3VER«R:~Sn  ^  »• 

lus^outint  *sl*.  r.,  %,  s> 

int«p»»*  *(10) 

omwntici’i  *<10,  10).  3(10) 

=  THIS  PORTION  DOES  FORWARD  DOSS 

so  10  »»£,  h 

»u»»0. 

so  20  j*l, »-l 

20  t uawtum*'*  (•  l  i  ) ,  j  )  *0  (•  ( J  >  1 

10  0 (» < l > > mO (* < l > > —sum 

C  THIS  PORTION  DOES  BODKUORD  P«SS 

b<s(n>  >»b(*(n> ) /«<»<n)  ,  n) 
do  30  i«rr-l,  1,-1 

*um>0. 

do  *0  j“i*l,n 

40  •  <•  <1 ) ,  J >*D <*  *  J }  * 

30  b(*<i>)“tb<»<l>  >  — *um) /•<»(!>, i> 

rot  urn 
•nd 


Y'*> 


128 


APPENDIX  B 

COMPUTER  PROGRAM  FOR  THE  BENDING  OF  A 
SINGLE  CANTILEVER  BEAM 
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c 

c 
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c 

c 

c 

c 


c 

c 

c 

c 


•  MAIN  PROGRAM  BEGINS 


I KTEGER  I EF 

REAL  DCADRE . AERR , RERR , ERROR . F , FI , NEW , CURR , CEE , EE I M 
REAE  BETA , E . S7I FF , P , V , D , PE , PB , ME , PEBAR , PBBAR , MBBAR 
REAE  A, B.C.DL. THETA, KSQ, CONST, EPS, XA , YA , NEWDEG , PEE 
REAL  LOLI K , OPEI K , MO 

EXTERNAL  FUNCTIONS  FOR  THE  INTEGRATION  SUBROUTINE  DCADRE. 

EXTERNAL  FEEN , XCOD , YCOD 

COMMON  BLOCKS  BETWEEN  FUCTIONS,  SUBROUTINES  AND  THE  MAIN 
PROGRAM  FOR  THE  DATA  TO  BE  COMMUNICATED. 

COMMON /AREA1 /BETA , L , STI FF , P , W , D , MBBAR 
COMMON  /AREA3/DL.C, THETA, KSQ, CONST 

DEFINE  ERROR  LIMITS  AND  TOLERANCES. 

EPS-1. OE-5 
AERR- 0.0 
RERR-1.0E-2 

READ  IN  GEOMETRY  AND  LOAD  ON  THE  ELEMENT. 

READ (1,*) BETA, L,STIFF,PLD,W,DBYL, MO 

OBTAIN  PRESSURE  LOAD  ECCENTRICITY  FROM  ECCENTRI Cl TY /LENGTH  RATIO. 
D-DBYL*L 

LOOP  TO  DETERMINE  TRAJECTORY  DUE  TO  VARYING  LEVELS  OF  PRESSURE. 

DO  81  K-1,7 
r-PLD*FLOAT(K) 

MB-P*D«-M0 
MBBAR-MB*L/STI FF 

INITIALLY  ASSUMED  EXTREME  ANGLES  FOR  THE  SECANT  METHOD. 

CURR-0./180.*3.14159 
NEW-175./180.*3 .14159 
CALL  COMP (CURR) 

UPLI M- ( THETA+CURR ) /2 . 

IF(1-KSQ*SIN(X)*SIN(X) )  >-0  SET  THE  UPPER  BOUND  SUCH  THAT 
THE  EQUALITY  HOLDS. 

UPTEMP-SQRT ( 1 . /KSQ ) 

I F ( UPTEMP . GT . 1 . )  THEN 
UPPER-1.57 

JLSE 

UPPER- AS I N ( UPTEMP ) 

ENDIF 

SET  THE  UPPER  INTEGRATION  LIMIT  TO  THE  LIMITING  UPPER  BOUND  IN 
EITHER  THE  POSITIVE  OR  NEGATIVE  RANGE. 

IF (UPLIM.GT. 0 . 0 )THEN 

JF  ( TTPT.TM .  GT .  TTPPF.P )  TTPT.T  M.nPPF.P 
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ELSE 

I F ( UPLI M . LT . “ UPPER ) UPLI M-  - UPPER 
ENDIF 

F-CONST*DCADRE(FLEN, THETA/2. , UPLIM.AERR, 

♦  RERR , ERROR , I ER ) ~1 . 

C  ITERATIVE  LOOP  TO  DETERMINE  THE  TIP  SLOPE  OF  THE  ELEMENT. 

C 

10  Fl-F 

OLD-CURR 
CURR-NEW 
CALL  COMP(CURR) 

UPLIM-(THETA*CURR)/2. 

C  IF(1“KSQ*SIN(X)*SIN(X)  >«  0  ,  SET  THE  UPPER  BOUND  SUCH  THAT 

C  THE  EQUALITY  HOLDS. 

C 

UPTEMP- SQRT  ( 1 .  /KSQ ) 

I F ( UPTEMP . GT . 1 . )  THEN 
UPPER-1.57 

ELSE 

UPPER-ASIN (UPTEMP ) 

ENDIF 

C  SET  UPPER  INTEGRATION  LIMIT  TO  THE  LIMITING  BOUND  IN  EITHER 

C  THE  POSITIVE  OR  THE  NEGATIVE  RANGE. 

IF (UPLIM.GT. 0 . 0 )THEN 

I F ( UPLIM . GT . UPPER ) UPLI M-UPPER 
ELSE 

IF (UPLIM.LT. “UPPER )UPLIM-~UPPER 
ENDIF 

F-CONST*DCADRE ( FLEN , THETA/2 . , UPLI M , AERR , RERR , 

♦  ERROR, IER)“1. 

C  IF  THE  DIFFERENCE  BETWEEN  THE  TWO  EXTREME  FUNCTIONS  IN  THE 

C  SECANT  METHOD  IS  <  EPS  THEN  WRITE  WARNING  MESSAGE  AND  QUIT 

C  TIP  ANGLE  COMPUTATION  LOOP.  USE  THE  LAST  COMPUTED  ANGLE 

C  ALPHA  AS  THE  TIP  ANGLE. 

C 

I F  ( ABS  ( F-Fl ) .  LT .  EPS )  THEN 
.  WRITE(3,23) 

23  FORMAT ( //5X , ’  F  AND  FI  DIFFER  BY  LESS  THAN  EPS’) 

GO  TO  20 
ENDIF 
C 

C  UPDATE  THE  ANGLE  ALPHA  USING  THE  SECANT  METHOD. 

C 

NEW-CURR-F* ( CURR-OLD ) / ( F-Fl ) 

C 

C  ERROR  CRITERION  TO  QUIT  TIP  ANGLE  COMPUTATION  LOOP. 

C 

IF  ( ABS  ( NEW-CURR ) . LT . EPS ) GO  TO  20 
WRITE (3 ,1 )NEW 

1  FORMAT (//’  NEW  ALPHA-’ , FI 0.4) 

GO  TO  10 

20  NEWDEG-NEW*180./3. 14159 

WRI TE ( 3 , 2 ) NEWDEG 

2  FORMAT (//’  ALPHA-’ ,F10.4,’  DEGREES'  ) 

CALL  COMP (NEW) 

UPLIM- ( THETA+NEW ) /2 . 

C 

C  IF(1-KSQ*SIN(X)*SIN(X) )  >-  0  ,  SET  THE  UPPER  BOUND  SUCH  THAT 

C  THE  EQUALITY  HOLDS. 

C 

UPTEMP- SORT ( 1 . /KSO ) 
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I F  ( UPTEMP . G7 . 1 .  1  THEN 
UPPER- 1.57 

ELSE 

UPPER- AS I N ( UPTEMP ) 

ENDIF 

C 

C  SET  UPPER  INTEGRATION  LIMIT  TO  THE  LIMITING  BOUND  IN  EITHER  THE 

C  POSITIVE  OR  NEGATIVE  RANGE. 

C 

IF (UPLIM.GT. 0 . 0 )THEN 

I F ( UPLI M . GT . UPPER ) UPLI M- UPPER 
ELSE 

IF (UPLIM.LT. -UPPER JUPLIM--UPPER 
ENDIF 

ANS-CONST*DCADRE ( FLEN , THETA/2 . , UPLI M , AERR , 

*  RERR , ERROR ,  I  ER ) 

WRI TE (3  7 )  ANS 

7  FORMAT (//’  ANSWER  OBTAINED  BY  INTEGRATION  IS  ’ , FI 0.4) 

LOLIM-THETA/2 . 

C 

C  COMPUTE  TIP  DEFLECTIONS  OF  THE  ELEMENT. 

C 

XA«CONST*DCADR£ ( XCOD , LOLIM , UPLI M , AERR , RERR , ERROR , I ER ) 
YA-CON5T*DCADR£ ( YCOD , LOLI M , UPL I M , AERR  f  RERR , ERROR , I ER ) 

WRITE( 3 ,12 ) 

WRITE (3,11) XA*L , YA*L 

12  FORMAT (//'  THE  NEW  COORDINATES  FOR  THE  TIP  IS  : '//) 

11  FORMAT ( '  XA- ' ,F1Q . 4 , '  YA- ’ ,F10 . 4//) 

WRITE ( 3,13)  . 

13  F0RMAT(//2X, 'NO. ’ ,92, 'THETA  (RADI ANS ) ’ , 82, ' XCOOD ' , 102 , 

«■  '  YCOOD'  //) 

C 

C  COMPUTE  DEFLECTIONS  OF  ELEMENT  NEARER  THE  TIP. 

C 

COUNT- 1. 

DO  21  J«l,10 

LLIM-NEW/ ( 2 . *COUNT ) +LOLIM 

XCOOD-CONST*DCADRE ( XCOD , LLI M , UPLI M , AERR , RERR , ERROR , I ER ) 
YCOOD-CONST*DCADRE ( YCOD , LLIK ,UPLIM , AERR , RERR , ERROR , I ER ) 
WRITE(3,9)J.LLIM,XC00D,YC00D 
9  FORMAT ( /I 3,5X,E20.6,5X,F10.4,5X,F10.4) 

COUNT-COUNT* 0 . 1 

21  CONTINUE 
C 

C  COMPUTE  DEFLECTIONS  OF  ELEMENT  NEARER  THE  BASE. 

C 

DO  22  J-3,6 

LLIM-NEW/ ( 2 . *FLOAT( J ) ) ■+LOLIM 

XCOOD-CONST*DCADR£ ( XCOD , LLI K , UPLI M , AERR , RERR , ERROR , I ER ) 

YCOOD- CONST*DCADRE ( YCOD , LLIM , UPLI K , AERR , RERR , ERROR , I ER ) 

M-J*8 

WRI TE ( 3 , 9 ) M , LLI M , XCOOD*  L , YCOOD*L 

22  CONTINUE 

81  CONTINUE 

STOP 

END 

C 

Q  ************************************************* 

C  *  FUNCTION  FLEN  TO  COMPUTE  THE  INTEGRAND  * 

Q  ************************************************* 

C 

FUNCTION  FLEN ( X ) 

REAL  2,KSQ 

COMMON  /AREA3/DL,C, THETA, KSQ, CONST 
DEFN-1 . -KSQ*SIN (X ) *SIN ( X ) 

ULLDEFN . LT . 0 . 0 ) DEFN- 1 . E- 4 


V- 

■I- 

> 

> 


* fc  . . ' a 


FLEN-l/SQRTIDEFN) 

RETURN 

END 


*  SUBROUTINE  COMP  TO  COMPUTE  THE 

*  INTERMEDIATE  VALUES  FOR  THE 

*  INTEGRATION 


SUBROUTINE  COMP (CURR  ) 

REAL  CURR , THETA , KSQ , CONST 
REAL  BETA.L, STIFF, P,W,D,MBBAR 
COMMON  /AREA1 /BETA , L . STIFF , P , V , D , MB BAR 
COMMON  /AREA3/DL.C, THETA, KSQ, CONST 
PL--P*COS  (  CURR )  «-W*COS  (  BETA ) 

PB-P*S I N ( CURR ) -W*  SI N ( BETA ) 

PLBAR-PL*L*L/STIFF 
PBBAR-PB*L*L/STI FF 
A-2.*PLBAR 
B-2 . *PBBAR 

C--2 . *PLBAR*COS ( CURR ) *2 . *PBBAR*SIN ( CURR ) *MBBAR*MBBAR 
DL-SQRT(A*A«-B*B) 

THETA-ATAN2 (B , A ) 

KSQ-2.*DL/(C+DL) 

CONST- 2 . /SQRT (C*DL) 

RETURN 

END 


**■**★*****•***★★*★•*★*•***★**  *★*'**'**  ******* 

*  FUNCTION  ZCOD  TO  COMPUTE  THE  * 

*  INTEGRAND  FOR  THE  COMPUTATION  * 

*  OF  THE  TIP  AND  INTERMEDIATE  * 

*  X-COORDINATE  * 

*****************  ************************ 

FUNCTION  ZCOD (2) 

REAL  Z 

COMMON  /AREA 3 /DL , C , THETA , KSQ , CONST 
ZCOD- COS ( 2 . *Z-THETA ) *FLEN ( Z ) 

RETURN 

END 

**************************************** 

*  FUNCTION  7C0D  TO  COMPUTE  THE  * 

*  INTEGRAND  FOR  THE  COMPUTATION  * 

*  OF  THE  TIP  AND  INTERMEDIATE  * 

*  7-COORDINATE  * 

******************************************* 

FUNCTION  7C0D(Z) 

REAL  Z 

COMMON  /AREA3/DL,C, THETA, KSQ, CONST 
7C0D- S I N ( 2 . *  Z-THETA ) *  FLEN ( Z ) 

RETURN 


t  / 


134 


C 

c  ..******».*t«*t....*.*..«,..,.**,*,* 

C  *  MAIN  PROGRAM  BEGINS  * 

c  **••**•**•**•*•••**«***•••**•«•«••,« 

C 

INTEGER  IER 

REAL  DCADRE , AERR , RERR , ERROR , F , FI , NEW , CURB , OLD , LLI K 
REAL  BETA , L , STI FF , P , W , D , PL , PB , MB , PLBAR , PBBAR , MBBAR 
REAL  A . B , C ,DL , THETA , KSQ , CONST , EPS , XA , YA , NEWDEG , PLD 
REAL  LOLIM , UPLIM ,M0 ,2(10) , Y ( 10 ) , ALPHA (10) , XXX (10) ,YYY(10) 

REAL  LEN(IO) ,EI (10) , I NTLO (10), SLFWT (10), DLRAT (10) 

REAL  EXTMOM(IO) 

C 

C  EXTERNAL  FUNCTIONS  FOR  THE  INTEGRATION  SUBROUTINE  DCADRE. 

C 

EXTERNAL  FLEN , XCOD , YCOD 
C 

C  COMMON  BLOCKS  BETWEEN  FUNCTIONS,  SUBROUTINES  AND  THE  MAIN  PROGRAM 

C  FOR  THE  DATA  TO  BE  COMMUNICATED. 

C 

COMMON /AREA1 /BETA , L, STIFF, PA, W,D, MBBAR 
COMMON  /AREA3/DL.C, THETA, XSQ, CONST 
C 

C  DEFINE  ERROR  LIMITS  AND  TOLERANCES. 

C 

EPS-1. OE-5 
AERR- 0.0 
RERR-1.0E-2 
DEL-0.2 
C 

C  INITIALLY  ASSUMED  GLOBAL  COORDINATES  OF  LOAD. 

C 

GLXP- 0.0 
GLYP-0.0 
C 

C  READ  IN  DATA  OF  LOAD  ,  LINK  JOINING  LIFTING  ELEMENTS  AND  END 

C  EFFECTOR  AND  THE  LAST  LIFTING  ELEMENT. 

C 

READ  ( 1 ,  *  )  WL ,  RL ,  XL  ,  YL ,  K ,  RB 
C 

C  READ  IN  TRANSFORMATION  MATRICES  FOR  EACH  OF  THE  LINKS  IN  THE 

C  LIFTING  ELEMENTS. 

C 

DO  82  K-2,6,2 

READ ( 1 , * ) X ( K ) , Y(K) , ALPHA ( K ) 

62  CONTINUE 

C 

C  READ  IN  GEOMETRY  AND  LOADS  ON  EACH  OF  THE  LIFTING  ELEMENTS ■ 

C 

DO  81  K-1,4 

^  READ ( 1 , * ) LEN ( K ) ,EI (K) , INTLD (K ) ,SLFWT(K) , DLRAT ( K ) , EXTMOM ( K ) 

C  ITERATIVE  LOOP  TO  DETERMINE  THE  GLOBAL  LOAD  POSITION  DUE  TO  A 

C  SET  OF  PRESSURES. 

C  SET  ANGLE  AT  WHICH  SHEAR  LOAD  AND  SELF  WEIGHT  ACTS  TO  180 

C  DEGREES  FOR  THE  FIRST  ELEMENT. 

C 

81  CONTINUE 

87  BETA-3.14159 

C 

C  LOOP  TO  DETERMINE  THE  LOCAL  (X,Y)  COORDINATES  FOR  THE  TIP  OF 

C  EACH  ELEMENT. 

C 


ALPSTTM-n.O 


STIFF-EI (K) 

PA-INTLD(K) 

SHEAR  LOAD  AND  SELF  WEIGHT  ACTING  ON  ELEMENT  K. 

W-WL-SLFWT(K) 

DBYL-DLRAT(K) 

MO-EXTMOM(K) 

D«DBYL*L 

MB-PA*D+M0 

MBBAR-MB*L/STIFF 

INITIALLY  ASSUMED  EXTREME  ANGLES  FOR  THE  SECANT  METHOD. 

CURR-0./180. *3.14159 
NEW-175 . /I 80. *3. 14159 
CALL  COMP(CURR) 

UPLIM-  ( THETA  "-CURR )  /2 . 

IF(1-KSQ*SIN(X)*SIN(X) )  >-  0  SET  THE  UPPER  BOUND  SUCH  THAT 
THE  EQUALITY  HOLDS. 

UPTEMP-SQRT ( 1 . /KSQ ) 

I F ( UPTEMP . GT . 1 . )  THEN 
UPPER-1.57 

ELSE 

UPPER- AS I N ( UPTEMP ) 

ENDIF 

SET  UPPER  INTEGRATION  LIMIT  TO  THE  LIMITING  UPPER  BOUND  IN 
EITHER  THE  POSITIVE  OR  NEGATIVE  RANGE. 

I F ( UPLI M.GT.0.0) THEN 

I F ( UPLIM . GT . UPPER ) UPLI M-UPPER 
ELSE 

I F ( UPLI M . LT . -UPPER ) UPLI M- -UPPER 
ENDIF 

LOLIM-THETA/2 . 

F-CONST*DCADRE ( FLEN , LOLIM , UPLIM , AERR , 

*  RERR, ERROR, IER)-1. 

ITERATIVE  LOOP  TO  DETERMINE  THE  TIP  ANGLE  ALPHA  IN  EACH  ELEMEN 
Fl-F 

OLD-CURR 
CURR-NEW 
CALL  COMP(CURR) 

UPLIM- ( THETA+CURR ) /2 . 

IF(1-KSQ*SIN(X)*SIN(X) )  >-  0  ,  SET  THE  UPPER  BOUND  SUCH  THAT 
THE  EQUALITY  HOLDS. 

UPTEMP-SQRT { 1 . /KSQ ) 

I F ( UPTEMP . GT . 1 . )  THEN 
UPPER-1.57 

ELSE 

UPPER-ASIN (UPTEMP) 

ENDIF 

SET  UPPER  INTEGRATION  LIMIT  TO  THE  LIMITING  BOUND  IN  EITHER 
THE  POSITIVE  OR  NEGATIVE  RANGE. 

IF (UPLIM. GT.0.0)THEN 

IF  ( ITPT.TM .  fiT  .  HPPF.R  )  T7PT.T  MmTTPPF.R 


nnno  ~J  nnnn  n  n  n  n  kimm  non  non  k>  nnnono 
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ELSE 

I F ( UPLI K.LT.- UPPER ) UPLI M- -UPPER 
ENDIF 

LOLIM-THETA/2. 

F -CONST*DCADRE ( FLEN , LOLIM  , UPLI M , AERR , RERR , 

-  ERROR, IER)-i. 

IF  THE  DIFFERENCE  BETWEEN  THE  TWO  EXTREME  FUNCTIONS  IN  THE 
SECANT  METHOD  IS  <  EPS  THEN  WRITE  WARNING  MESSAGE  AND  QUIT 
TIP  ANGLE  COMPUTATION  LOOP.  USE  THE  LAST  COMPUTED  ANGLE 
ALPHA  AS  THE  TIP  ANGLE. 

IF ( ABS (F-Fl ) .LT.EPS )THEN 
WRITE( 3,23) 

FORMAT ( //52, 'F  AND  FI  DIFFER  BY  LESS  THAN  EPS’) 

GO  TO  20 
ENDIF 

UPDATE  THE  ANGLE  ALPHA  USING  THE  SECANT  METHOD. 

NEW-CURR-F* ( CURR-OLD ) / ( F-Fl ) 

ERROR  CRITERION  TO  QUIT  TIP  ANGLE  COMPUTATION  LOOP. 

IF  (ABS (NEW-CURR) .LT.EPS) GO  TO  20 
WRI TE (3,1) NEW 

FORMAT (//'  NEW  ALPHA-* , FI 0.4) 

GO  TO  10 

NEWDEG-NEW-180 ./3 .14159 
WRI TE (3,2) NEWDEG 

FORMAT ( // '  ALPHA-’ ,F10.4, *  DEGREES’  ) 

CALL  COMP (NEW) 

UPLIM- ( THETA+NEW ) /2 . 

IF(1-KSQ*SIN(X)*SIN(S) )  >•  0  ,  SET  THE  UPPER  BOUND  SUCH  THAT 
THE  EQUALITY  HOLDS. 

UPTEMP-SQRT ( 1 . /KSQ ) 

I F ( UPTEMP . GT . 1 . )THEN 
UPPER-1.57 

ELSE 

UPPER-AS I N ( UPTEMP ) 

ENDIF 

SET  UPPER  INTEGRATION  LIMIT  TO  THE  LIMITING  BOUND  IN  EITHER  THE 
POSITIVE  OR  NEGATIVE  RANGE. 

IF (UPLIM. GT.0.0)THEN 

I F ( UPLI M . GT . UP  PER ) UPL I M- UPPER 
ELSE 

I F ( UPLIM.LT . -UPPER ) UPLI M- -UPPER 
ENDIF 

LOLIM-THETA/2. 

ANS-CONST*DCADRE ( FLEN , LOLIM , UPLIM , AERR , 

*  RERR, ERROR, I ER) 

WRITE (3, 7  JANS 

FORMAT (//'  ANSWER  OBTAINED  BY  INTEGRATION  IS  ’,F10.4) 
LOLIM-THETA/2. 

XA - CON ST*DCADRE ( XCOD , LOLI M , UPLIM , AERR , RERR , ERROR , I ER ) 
YA-CONST*DCADRE ( YCOD , LOLIM , UPLI M , AERR , RERR , ERROR , I ER ) 

STORE  TIP  DEFLECTIONS  AND  CLOCKWISE  ANGLE  ROTATION  IN 
ARRAYS  X,Y  AND  ALPHA. 

X(2*K-1)-XA*L 
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ALPHA  (2*K-1)  —NEW 


TOTAL  CHANGE  IN  ANGLE  FROM  THE  SUPPORTING  BASE. 


ALPSUM-ALPSUM*NEW 

WRITE  THE  TIP  DEFLECTIONS  OF  THE  ELEMENT  IN  THE  LOCAL  COORDINATE 
SYSTEM. 


WRITE ( 3,21)X(2*K-1) , Y ( 2*K-1 ) 

FORMAT ( //5X, ’ZTIP-' ,F7 . 3 , 5X, * YTIP- * ,F7 . 3 ) 

UPDATE  THE  ANGLE  AT  WHICH  THE  DEAD  LOAD  ACTS. 

I F ( K . NE . 4 ) BETA- BETA-NEW 
CONTINUE 


CALL  TO  TRANSFORM  TIP  DEFLECTIONS  IN  LOCAL  COORDINATE  SYSTEM 
TO  GLOBAL  COORDINATE  SYSTEM. 

CALL  TRANS (X, Y, ALPHA, XXX, YYY, 4 5 

WRITE  GLOBAL  TIP  DEFLECTIONS  OF  ELEMENTS. 

DO  83  K-1,4 

WRITE(3,115)K,XXX(K) ,K,TYY(K) 

CONTINUE 

FORMAT ( //5X,  ’GLTPXELC ,12 , ' )- ' ,F7 . 3 , 5Z, ' GLTPYEL ( ',12, ’ , 

*  F7.3) 


COMPUTE  COORDINATES  OF  LOAD. 


i 


3 


A 


RTPL-2 . •RB+H+RL 
PS I -ATAN2 { XL . RL+YL ) 

GLXC-XXX ( 4  )  -*-RTPL*SIN  { ALPSUM+P5I ) 
GLYC-YTY  { 4 ) -RTPL*COS ( ALPSUM+PSI ) 


CRITERION  TO  QUIT  GLOBAL  LOAD  POSITION  ITERATION  DUE  TO 
A  SET  OF  PRESSURES. 


IF ( ( ABS ( GLXC-GLXP } .LE.DEL) .AND. ( ABS (GLYC-GLYP) .LE.DEL ) )GO  TO  86 
UPDATE  EXTERNAL  MOMENT  IN  EACH  ELEMENT  DUE  TO  LOAD. 

DO  84  1-1,4 

EXTMOM ( I ) - ( SLFWT ( I ) +WL ) * ( YYY ( I ) -GLYF ) 

CONTINUE 

UPDATE  GLOBAL  LOAD  POSITION. 

GLXP-GLXC 
GLYP-GLYC 
GO  TO  87 

WRITE  FINAL  GLOBAL  LOAD  POSITION  DUE  TO  A  SET  OF  PRESSURES. 
WRITE ( 3,22 )GLXC ,GLYC 

FORMAT ( //5X, 'GL.  XCOOD.  LD . - ’ , F10 . 4 , 5X , ' GL.  YCOOD.  LD.-’,F10.4) 

STOP 

END 


*********************♦*#**#***##***#********«**  +  ***'** 

*  FUNCTION  FLEN  TO  COMPUTE  THE  INTEGRAND  * 

W************************  *****************  *********** 


JUNCTION  FLEN(X) 
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REAL  X,  KSQ 

COMMON  /AREA3/DL,C, THETA, KSQ, CONST 
DEFN-1 . “KSQ*SIN (X ) *SIN (2 ) 

IF (DEFN .LE. 0 . 0 )DEFN-1 . OE-4 
FLEN«1./SQRT(DEFN) 

RETURN 

END 


SUBROUTINE  COMP  TO  COMPUTE  THE 
INTERMEDIATE  VALUES  FOR  THE 
INTEGRATION 


SUBROUTINE  COMP(CURR) 

REAL  CURR, THETA, KSQ, CONST 
REAL  BETA, L, STIFF, P,W,D,MBBAR 
COMMON  /AREA1/BETA , L , STI FF , P , K , D , MBBAR 
COMMON  /AREA3/DL , C , THETA , KSQ , CONST 
PL— P*COS  ( CURR )  ♦W»COS  ( BETA ) 

PB-P*SIN { CURR ) -W*SIN ( BETA ) 

PLBAR-  PL  *L*L/STI FF 
PBBAR«PB*L*L/STI FF 
A-2.*PLBAR 
B-2.*PBBAR 

C— 2 . * PLBAR* COS ( CURR) +2 . *PBBAR*SIN ( CURR ) +MBBAR*M3BAR 
DL»SQRT(A*A+B*B) 

TH£TA-ATAN2(B,A) 

KSQ-2.*DL/(C+DL) 

CONST- 2 . /SQRT ( C+DL ) 

RETURN 

END 


FUNCTION  XCOD  TO  COMPUTE  THE 
INTEGRAND  FOR  THE  COMPUTATION 
OF  THE  TIP  X-COORDINATE 


C***t*t****** 


FUNCTION  XCOD(X) 

REAL  X 

COMMON  /AREA3/DL,C, THETA, KSQ, CONST 
XCOD-COS ( 2 . *X-THETA ) *FLEN ( X ) 

RETURN 

END 


**************************************** 

*  FUNCTION  YCOD  TO  COMPUTE  THE  * 

*  INTEGRAND  FOR  THE  COMPUTATION  * 

*  OF  THE  TIP  Y-COORDINATE  * 

**************************************** 


FUNCTION  YCOD(X) 

REAL  X 

COMMON  /AREA3/DL,C, THETA, KSQ, CONST 
YCOD-SIN ( 2 , *X-THETA ) *FLEN (X ) 

RETURN 

END 

**■******«*************•»****«*«"********** 

*  SUBROUTINE  TRANS  TO  TRANSFORM  * 

*  THE  LOCAL  COORDINATES  TO  THE  * 

*  GLOBAL  COORDINATES  * 
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SUBROUTINE  TRANS  (2,  Y , ALPHA , 222 , YYY , I ELEM ) 

REAL  2(10)  ,Y(10)  ,  ALPHA  (10)  ,  TRAN  (3, 3)  ,  TEMP  (3, 3)  ,  TRAN  2(3,3) 
REAL  222(10) , YYY (10) 

C 

C  INITIALIZE  TRANSFORMATION  MATRI2  TO  AN  IDENTITY  MATRI2. 

C 

DO  11  1-1,3 
DO  12  J-l,  3 

IF  (I.EQ.J)  THEN 
TRAN (I , J)-l. 

ELSE 

TRAN (I , J ) -0  . 

ENDIF 

12  CONTINUE 

11  CONTINUE 

C 

C  GLOBAL  AND  LOCAL  COORDINATES  FOR  FIRST  ELEMENT  COINCIDES. 

C 

222(1 ) -X( 1 ) 

YYY(1)-Y(1) 

DO  15  J-2 , I ELEM 
DO  10  I-2*J-3 , 2* (J-l ) 

CALL  UTOI ( X  ( I  )  ,  Y  ( I  )  , ALPHA ( I ) , TEMP ) 

CALL  TMULT< TRAN, TEMP, TRAN 2 ) 

C 

C  UPDATE  TRANSFORMATION  MATRIX  TILL  ELEMENT (J-l)  AND 

C  THE  RIGID  LINK  AFTER  IT. 

C 

DO  51  11-1,3 
DO  50  JJ-1,3 

TRAN ( 1 1 ,  J  J ) -TRAN2 ( I I ,JJ) 

50  CONTINUE 

51  CONTINUE 

10  CONTINUE 

C 

C  COMPUTE  GLOBAL  TIP  COORDINATES  FOR  ELEMENT  J. 

C 

222 ( J ) -TRAN ( 1 , 1 ) *2 ( 2* J-l ) +TRAN ( 1 , 2 ) *Y ( 2* J-l ) -TRAN ( 1 , 3 ) 

YYY ( J ) -TRAN ( 2 , 1 ) *2 ( 2* J-l ) -TRAN ( 2 , 2 ) *Y ( 2* J-l ) -TRAN ( 2 , 3 ) 

15  CONTINUE 

RETURN 
END 
C 

c  **************************************** 

C  *  SUBROUTINE  UTOI  TO  FORM  THE  * 

C  *  TRANSFORMATION  MATRIX  BETWEEN  * 

C  *  FRAME  (I)  AND  FRAME  (1-1)  * 

c  **************************************** 

C 

SUBROUTINE  UTOI ( 2 , Y , ALPHA , I  FORM) 

REAL  2 , Y , ALPHA , I FORM (3,3) 

I FORM ( 1 , 1 ) -COS ( ALPHA ) 

IFORM( 2,1) --SIN (ALPHA) 

I FORM ( 3 , 1 )-0 . 0 
IFORM( 1,2) -SIN (ALPHA) 

I FORM ( 2 , 2 ) -COS (ALPHA ) 

I FORM ( 3 , 2 )-0 . 0 
I FORM ( 1 , 3 )-X 
I FORM ( 2 , 3 ) «  Y 
I FORM ( 3 , 3 ) -1 . 0 
RETURN 
END 
C 

Q  **************************************** 

C  *  SUBROUTINE  TMULT  TO  MULTIPLY  * 

C  *  TWO  MATRICES  * 


SUBROUTINE  TMULT ( BREL A , CRELB , OREL A ) 

REAL  BRELA (3,3), CRELB (3,3), CRELA (3,3) 

DO  10  1-1,3 
DO  20  J-1,3 

CRELA C , J ) -BRELA C , 1 ) *CRELB ( 1 , J > -BRELA ( I ,2 '* 
CRELB ( 2 , J ) -BRELA (1,3)*  CRELB ( 3 , J ) 

CONTINUE 
CONTINUE 
RETURN 


